CN113093242B - 一种基于球谐展开的gnss单点定位方法 - Google Patents

一种基于球谐展开的gnss单点定位方法 Download PDF

Info

Publication number
CN113093242B
CN113093242B CN202110283670.9A CN202110283670A CN113093242B CN 113093242 B CN113093242 B CN 113093242B CN 202110283670 A CN202110283670 A CN 202110283670A CN 113093242 B CN113093242 B CN 113093242B
Authority
CN
China
Prior art keywords
observation
value
formula
satellite
iteration
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202110283670.9A
Other languages
English (en)
Other versions
CN113093242A (zh
Inventor
郭金运
郭恒洋
杨洲铭
邢云鹏
刘新
孔巧丽
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Shandong University of Science and Technology
Original Assignee
Shandong University of Science and Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Shandong University of Science and Technology filed Critical Shandong University of Science and Technology
Priority to CN202110283670.9A priority Critical patent/CN113093242B/zh
Priority to CN202210137675.5A priority patent/CN114518586B/zh
Publication of CN113093242A publication Critical patent/CN113093242A/zh
Priority to PCT/CN2021/123845 priority patent/WO2022048694A1/zh
Application granted granted Critical
Publication of CN113093242B publication Critical patent/CN113093242B/zh
Priority to ZA2022/03722A priority patent/ZA202203722B/en
Priority to US17/730,567 priority patent/US20220299652A1/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/01—Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/03—Cooperating elements; Interaction or communication between different cooperating elements or between cooperating elements and receivers
    • G01S19/07—Cooperating elements; Interaction or communication between different cooperating elements or between cooperating elements and receivers providing data for correcting measured positioning data, e.g. DGPS [differential GPS] or ionosphere corrections
    • G01S19/072—Ionosphere corrections
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/38—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system
    • G01S19/39—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system the satellite radio beacon positioning system transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/42—Determining position
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/01—Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/03—Cooperating elements; Interaction or communication between different cooperating elements or between cooperating elements and receivers
    • G01S19/07—Cooperating elements; Interaction or communication between different cooperating elements or between cooperating elements and receivers providing data for correcting measured positioning data, e.g. DGPS [differential GPS] or ionosphere corrections
    • G01S19/073—Cooperating elements; Interaction or communication between different cooperating elements or between cooperating elements and receivers providing data for correcting measured positioning data, e.g. DGPS [differential GPS] or ionosphere corrections involving a network of fixed stations
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/01—Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/03—Cooperating elements; Interaction or communication between different cooperating elements or between cooperating elements and receivers
    • G01S19/08—Cooperating elements; Interaction or communication between different cooperating elements or between cooperating elements and receivers providing integrity information, e.g. health of satellites or quality of ephemeris data
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/01—Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/13—Receivers
    • G01S19/24—Acquisition or tracking or demodulation of signals transmitted by the system
    • G01S19/25—Acquisition or tracking or demodulation of signals transmitted by the system involving aiding data received from a cooperating element, e.g. assisted GPS
    • G01S19/258—Acquisition or tracking or demodulation of signals transmitted by the system involving aiding data received from a cooperating element, e.g. assisted GPS relating to the satellite constellation, e.g. almanac, ephemeris data, lists of satellites in view
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/38—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system
    • G01S19/39—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system the satellite radio beacon positioning system transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/42—Determining position
    • G01S19/43—Determining position using carrier phase measurements, e.g. kinematic positioning; using long or short baseline interferometry
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/38—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system
    • G01S19/39—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system the satellite radio beacon positioning system transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/42—Determining position
    • G01S19/43—Determining position using carrier phase measurements, e.g. kinematic positioning; using long or short baseline interferometry
    • G01S19/44—Carrier phase ambiguity resolution; Floating ambiguity; LAMBDA [Least-squares AMBiguity Decorrelation Adjustment] method

Landscapes

  • Engineering & Computer Science (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Computer Security & Cryptography (AREA)
  • Position Fixing By Use Of Radio Waves (AREA)

Abstract

本发明属于卫星导航定位技术领域，具体公开了一种基于球谐展开的GNSS单点定位方法，该方法采用球谐展开描述与测站和卫星之间的高度角和方位角有关的误差项，包括对流层延迟误差以及电离层延迟误差，使用GNSS观测数据和卫星星历数据即可实现测站空间位置的计算。相比于现有的利用经验模型进行对流层延迟改正方法，本发明方法能够高效、快速地获得测点的位置信息，具有实际操作简单、数据处理便捷以及计算效率高等优势。

Description

一种基于球谐展开的GNSS单点定位方法
技术领域
本发明属于卫星导航定位技术领域，涉及一种基于球谐展开的GNSS单点定位方法。
背景技术
GNSS单点定位方法，是指根据卫星星历给出的观测瞬间卫星在空间的位置和卫星钟差，以及由一台GNSS接收机所测定的从卫星到GNSS接收机之间的距离，通过距离交会的方法，独立测定该GNSS接收机在地球坐标系中的三维坐标的定位方法。
利用广播星历提供的卫星位置和卫星钟差，以及伪距观测值进行的定位方式，称为标准单点定位。由于广播星历和伪距观测值精度的限制，且其数学模式中各项改正误差的影响又无法完全消除，因而，标准单点定位的精度通常只能达到米级。因此，标准单点定位主要用于导航定位及资源调查、勘测等低精度定位领域。利用精密星历和精密卫星钟差数据，由相位观测数据进行的定位方式，称为精密单点定位，精密单点定位的精度可以达到厘米量级。
在影响单点定位结果的三类误差源中，与信号传播有关的误差(主要是指电离层延迟和对流层延迟误差)无法很好地被模型化，因此，对定位的影响较大。
对于电离层延迟，单频接收机需要通过模型加以改正；双频接收机则可以通过线性组合法消除电离层延迟的影响；对于对流层延迟，由于对流层的性质为非弥散性，故不能用双频观测数据线性组合的方法消除对流层延迟误差，其对定位精度的影响不可以忽略。
对流层延迟分为干延迟和湿延迟两部分。在通常的定位解算策略中，天顶方向的干延迟由对流层延迟模型计算得出，通过映射函数将其映射到传播路径；而天顶方向的湿延迟则作为未知参数求解，并通过映射函数映射到传播路径。GNSS测量中对流层延迟模型或方法有萨斯塔莫宁模型、UNB3模型、气象元素法等。
以上三种对流层延迟改正方法存在如下问题：萨斯塔莫宁模型的精度相对较高，但是需要用到实测的气象数据，气象数据的获取问题，限制了萨斯塔莫宁模型的应用；UNB3模型不需要实测气象数据注入，但是UNB3模型的精度有限；气象元素法需要用到外部气象设备采集数据，不同设备的精度存在差异且设备笨重、昂贵，增加了外业采集数据的工作量；同时气象元素法也降低了单点定位的工作效率，很难为全天候的单点定位提供所需数据。
以上三种对流层延迟改正方法构建多种模型，以适应不同条件下改正对流层延迟的需要，考虑到实用性和普适性会使得测量外业和内业工作量增大，数据处理效率降低。此外，对于单点定位来说，对流层的延迟误差校正多采用经验模型，此类经验模型的参数估计也多采用经验值，因而，难以符合实际单点定位中的参数估计。
发明内容
本发明的目的之一在于提出一种基于球谐展开的GNSS标准单点定位方法，该方法基于球谐展开表示与测站和卫星间的方位角和高度角有关的误差项，以便能够快速、高效、准确地进行标准单点定位。本发明为了实现上述目的，采用如下技术方案：
一种基于球谐展开的GNSS标准单点定位方法，包括如下步骤：
I.1.建立利用伪距观测值进行的标准单点定位观测方程，如公式(1)所示：
Figure BDA0002979540070000021
式中，i表示第i观测历元，j表示卫星编号；ρ表示伪距观测值，
Figure BDA0002979540070000022
表示卫星j在第i观测历元的伪距观测值；
Figure BDA0002979540070000023
表示接收机的位置与卫星j在第i观测历元的位置之间的几何距离；c表示光速；tr表示测站接收机钟差，(tr)i表示第i观测历元测站的接收机钟差；ts表示卫星钟差，
Figure BDA0002979540070000024
表示卫星j在第i观测历元的钟差；
Figure BDA0002979540070000025
表示卫星j在第i观测历元的信号传播路径上的对流层延迟误差；
Figure BDA0002979540070000026
表示卫星j在第i观测历元的信号传播路径上的电离层延迟误差；ερ表示伪距观测数据残差；
定义卫星j在第i观测历元观测瞬间的空间三维坐标为
Figure BDA0002979540070000027
测站的近似坐标为(X0,Y0,Z0)，则卫星到测站近似位置的几何距离
Figure BDA0002979540070000028
表示为：
Figure BDA0002979540070000029
I.2.基于球谐展开表示与测站和卫星之间的高度角和方位角有关的误差项，误差项包括对流层延迟误差以及电离层延迟误差；基于球谐展开的标准单点定位观测方程表示为：
Figure BDA00029795400700000210
式中，n为球谐展开的阶数，m为球谐展开的次数，Nmax为球谐展开的最大阶数；
Figure BDA00029795400700000211
表示n阶m次的缔合勒让德多项式；Cnm和Snm分别表示n阶m次的球谐展开的系数，Cnm和Snm为球谐展开的待求参数；
Figure BDA00029795400700000212
分别表示测站与卫星j在第i观测历元之间的高度角和方位角；为了方便表示球谐展开，令：
Figure BDA0002979540070000031
则公式(2)简化为：
Figure BDA0002979540070000032
由简化后的标准单点定位观测方程(3)，得到标准单点定位误差方程(4)，公式如下：
Figure BDA0002979540070000033
其中，vρ表示伪距观测数据的改正值；
Figure BDA0002979540070000034
表示卫星j在第i观测历元的伪距观测数据的改正值；
I.3.基于公式(4)得到标准单点定位误差方程的线性化表达式，并进行简化处理，然后对简化后的标准单点定位误差方程的线性化表达式进行滑动解算；
I.3.1.对标准单点定位误差方程(4)在测站的近似坐标(X0,Y0,Z0)处进行泰勒级数展开并保留一阶项，则标准单点定位误差方程的线性化表达式如公式(5)所示；
Figure BDA0002979540070000035
公式(5)中，令：
Figure BDA0002979540070000036
Figure BDA0002979540070000037
Figure BDA0002979540070000038
其中，
Figure BDA0002979540070000039
表示测站近似坐标与卫星j在第i观测历元坐标计算出的待求参数dX的系数，待求参数dX是测站近似坐标X0的改正数；
Figure BDA00029795400700000310
表示测站近似坐标与卫星j在第i观测历元坐标计算出的待求参数dY的系数，待求参数dY是测站近似坐标Y0的改正数；
Figure BDA00029795400700000311
表示测站近似坐标与卫星j在第i观测历元坐标计算出的待求参数dZ的系数，待求参数dZ是测站近似坐标Z0的改正数；
Figure BDA00029795400700000312
表示由测站近似坐标与卫星j在第i观测历元坐标计算出的球谐展开的待求参数Cnm的系数；
Figure BDA00029795400700000313
表示由测站近似坐标与卫星j在第i观测历元坐标计算出的球谐展开的待求参数Snm的系数；d(tr)i表示在第i观测历元测站接收机钟差的待求参数；
则简化后的标准单点定位误差方程的线性化表达式，如公式(6)所示；
Figure BDA0002979540070000041
公式(6)中，i＝1,2,…,epoch，epoch表示最大观测历元数；
Figure BDA0002979540070000042
表示由卫星j在第i观测历元伪距观测值、卫星j在第i观测历元的钟差以及测站近似坐标与卫星j在第i观测历元之间的几何距离计算出的常数项，其中，
Figure BDA0002979540070000043
设定滑动窗口选择i个历元，每个历元最多能够观测到j颗卫星，则在这个滑动窗口下所有可视卫星的伪距观测值组成imax个方程，其中，imax＝i×j；
定义标准单点定位误差方程的系数矩阵、伪距观测值的改正数向量矩阵、常数项矩阵以及未知参数矩阵分别为矩阵B、V、L、X，则它们的表达式分别如下：
Figure BDA0002979540070000044
Figure BDA0002979540070000045
Figure BDA0002979540070000046
X＝[dX dY dZ (dtr)1 … (dtr)i C00 … CNmaxNmax SNmaxNmax]T；
则标准单点定位误差方程的线性化表达式(6)用一个矩阵形式表达，如公式(7)所示：
V＝BX-L (7)
未知参数X的最小二乘估值为：X＝(BTPB)-1BTPL (8)
公式(8)中，P为单位权矩阵；
I.3.2.判断公式(7)中系数矩阵B是否为病态，经过判断，若系数矩阵B是病态，则执行步骤I.3.3；否则，系数矩阵B是良态，执行步骤I.3.4；
I.3.3.首先利用截断奇异值分解的方法，为最小二乘谱修正迭代提供待求解的未知参数X的初值，然后基于最小二乘谱修正迭代计算未知参数X的值；
设系数矩阵B∈Rn×m，Rn×m表示n行m列的实数阵；则系数矩阵B的奇异值分解为：
B＝USVT (9)
公式(9)中，U∈n×n，V∈m×m，U、V均为正交矩阵，S∈n×m为一对角矩阵；
系数矩阵B∈Rn×m的截断奇异值矩阵Bk定义为：
Figure BDA0002979540070000051
公式(10)中，矩阵S中最小的(r-k)个非奇异值用零代替，即被截断，且k≤r；r表示系数矩阵B的秩，k表示矩阵S中保留的奇异值数量；
ui表示矩阵U对应的向量，vi表示矩阵V对应的向量，σi表示矩阵S中保留的奇异值；
计算矩阵S中奇异值的均值，将该均值作为截断奇异值的阈值，其中，矩阵S中奇异值大于截断奇异值的阈值的保留，小于截断奇异值的阈值的进行零化处理；
则公式(7)的截断奇异值解
Figure BDA0002979540070000052
为：
Figure BDA0002979540070000053
公式(11)中，
Figure BDA0002979540070000054
根据最小二乘原理VTPV＝min，公式(8)写成以下形式：
(BTPB)X＝(BTPL) (12)
将公式(12)的等式左端和右端同时加未知参数KX，并化简得到公式(13)：
(BTPB+KI)X＝BTPL+KX (13)
由公式(13)整理得，最小二乘谱修正迭代公式为：
X＝(BTPB+KI)-1(BTPL+KX) (14)
公式(14)中，K为任意实数；
基于最小二乘谱修正迭代计算未知参数X的解，转到步骤I.3.5；
I.3.4.利用最小二乘法求解未知参数X的解，转到步骤I.3.5；
I.3.5.根据计算得到的未知参数X的解，得到未知参数X中包含的参数值，即测站的位置参数(dX,dY,dZ)、接收机钟差dtr以及球谐展开的系数Cnm和Snm；
将未知参数X中包含的位置参数(dX,dY,dZ)，分别加到测站的近似坐标(X0,Y0,Z0)上，即标准单点定位求出的地心地固坐标系下的测站坐标(X,Y,Z)；
I.3.6.将测站坐标由地心地固坐标系转换到站心坐标系，实现GNSS标准单点定位。
本发明的目的之二在于提出一种基于球谐展开的GNSS精密单点定位方法，该方法基于球谐展开表示与测站和卫星间的方位角和高度角有关的误差项，以便能够快速、高效、准确地进行精密单点定位。本发明为了实现上述目的，采用如下技术方案：
一种基于球谐展开的GNSS精密单点定位方法，包括如下步骤：
II.1.建立利用载波相位观测值进行的精密单点定位观测方程，如公式(1)所示：
Figure BDA0002979540070000061
公式(1)中，i表示第i观测历元；j表示卫星编号；
Figure BDA0002979540070000062
表示载波相位观测值，即等效距离，
Figure BDA0002979540070000063
表示卫星j在第i观测历元的载波相位观测值，即等效距离；
Figure BDA0002979540070000064
表示接收机的位置与卫星j在第i观测历元的位置之间的几何距离；tr表示测站接收机钟差，(tr)i表示第i观测历元测站接收机钟差；ts表示卫星钟差，
Figure BDA0002979540070000065
表示卫星j在第i观测历元的钟差；
Figure BDA0002979540070000066
表示卫星j在第i观测历元的信号传播路径上的对流层延迟误差；
Figure BDA0002979540070000067
表示卫星j在第i观测历元的信号传播路径上的电离层延迟误差；N表示卫星相位观测数据的整周模糊度参数；
Figure BDA0002979540070000068
表示卫星j在第i观测历元的相位观测数据整周模糊度参数；c表示光速；λ表示载波的波长；
Figure BDA00029795400700000616
表示载波相位观测数据残差；
定义卫星j在第i历元观测瞬间的空间三维坐标为
Figure BDA0002979540070000069
测站的近似坐标为(X0,Y0,Z0)，则卫星到测站近似位置的几何距离
Figure BDA00029795400700000610
表示为：
Figure BDA00029795400700000611
II.2.基于球谐展开表示与测站和卫星间的高度角和方位角有关的误差项，误差项包括对流层延迟误差以及电离层延迟误差；基于球谐展开的精密单点定位观测方程表示为：
Figure BDA00029795400700000612
式(2)中，n为球谐展开的阶数，m为球谐展开的次数，Nmax为球谐展开的最大阶数；
Figure BDA00029795400700000613
表示n阶m次的缔合勒让德多项式；Cnm和Snm表示n阶m次的球谐展开的系数，Cnm和Snm为球谐展开待求参数；
Figure BDA00029795400700000614
和
Figure BDA00029795400700000615
分别表示测站与卫星j在第i观测历元之间的高度角和方位角；为了方便表示球谐展开，令：
Figure BDA0002979540070000071
则基于球谐展开的精密单点定位观测方程(2)简化为公式(3)；
Figure BDA0002979540070000072
由公式(3)中精密单点定位观测方程得到精密单点定位误差方程，如公式(4)所示；
Figure BDA0002979540070000073
其中，
Figure BDA0002979540070000074
表示载波相位观测数据的改正值；
Figure BDA0002979540070000075
表示卫星j在第i观测历元的载波相位观测数据的改正值；
II.3.基于公式(4)得到精密单点定位误差方程的线性化表达式，并进行简化处理，然后对简化后的精密单点定位误差方程的线性化表达式进行滑动解算；
I.3.1.对精密单点定位误差方程(4)在测站的近似坐标(X0,Y0,Z0)处进行泰勒级数展开并保留一阶项，则精密单点定位误差方程的线性化表达式如公式(5)所示；
Figure BDA0002979540070000076
公式(5)中，
Figure BDA0002979540070000077
表示卫星j在第i历元的相位观测数据模糊度初值；
Figure BDA0002979540070000078
表示整周模糊度参数的改正数；d(tr)i表示在第i观测历元测站接收机钟差的待求参数；
公式(5)中，令：
Figure BDA0002979540070000079
Figure BDA00029795400700000710
Figure BDA00029795400700000711
Figure BDA00029795400700000712
为测站近似坐标与卫星j在第i观测历元坐标计算出的待求参数dX的系数，待求参数dX是测站近似坐标X0的改正数；
Figure BDA00029795400700000713
为测站近似坐标与卫星j在第i观测历元坐标计算出的待求参数dY的系数，待求参数dY是测站近似坐标Y0的改正数；
Figure BDA00029795400700000714
为测站近似坐标与卫星j在第i观测历元坐标计算出的待求参数dZ的系数，待求参数dZ是测站近似坐标Z0的改正数；Anm表示球谐展开待求参数Cnm的系数，Bnm表示球谐展开待求参数Snm的系数；
Figure BDA0002979540070000081
表示由测站近似坐标与卫星j在第i观测历元坐标计算出的球谐展开的待求参数Cnm的系数，
Figure BDA0002979540070000082
表示由测站近似坐标与卫星j在第i观测历元坐标计算出的球谐展开的待求参数Snm的系数；
II.3.1.首先给出简化后的精密单点定位误差方程的线性化表达式，如公式(6)所示：
Figure BDA0002979540070000083
公式(6)中，i＝1,2,,…,epoch，epoch表示最大观测历元；
Figure BDA0002979540070000084
表示由卫星j在第i观测历元载波相位观测值、卫星j在第i观测历元的钟差、卫星j在第i历元的相位观测数据整周模糊度初始值以及测站近似坐标与卫星j在第i观测历元之间的几何距离计算出的常数项，
Figure BDA0002979540070000085
设定滑动窗口选择i个历元，每个历元最多能够观测到j颗卫星，则在这个滑动窗口下所有可视卫星的载波相位观测值组成imax个方程，其中，imax＝i×j；
定义精密单点定位误差方程的系数矩阵、载波相位观测值的改正数向量矩阵、常数项矩阵以及未知参数矩阵分别为矩阵B、V、L、X，则它们的表达式分别如下：
Figure BDA0002979540070000086
Figure BDA0002979540070000087
Figure BDA0002979540070000088
Figure BDA0002979540070000089
Figure BDA00029795400700000810
表示第1颗卫星在第1个观测历元的整周模糊度参数，
Figure BDA00029795400700000811
表示第j颗卫星在第i个观测历元的整周模糊度参数；
则精密单点定位误差方程的线性化表达式(6)用一个矩阵形式表达，如公式(7)所示；
V＝BX-L (7)
未知参数X的最小二乘估值为：X＝(BTPB)-1BTPL (8)
公式(8)中，P为单位权矩阵；
II.3.2.判断公式(7)的系数矩阵B是否为病态，若系数矩阵B是病态，则执行步骤II.3.3；否则，若系数矩阵B为良态，执行步骤II.3.4；
II.3.3.首先利用截断奇异值分解的方法，为最小二乘谱修正迭代提供待求解的未知参数X的初值，然后基于最小二乘谱修正迭代计算未知参数X的值；
设系数矩阵B∈Rn×m，Rn×m表示n行m列的实数阵；则系数矩阵B的奇异值分解为：
B＝USVT (9)
公式(9)中，U∈n×n，V∈m×m，U、V均为正交矩阵，S∈n×m为一对角矩阵；
系数矩阵B∈Rn×m的截断奇异值矩阵Bk定义为：
Figure BDA0002979540070000091
公式(10)中，矩阵S中最小的(r-k)个非奇异值用零代替，即被截断，且k≤r；r表示系数矩阵B的秩，k表示矩阵S中保留的奇异值数量；
ui表示矩阵U对应的向量，vi表示矩阵V对应的向量，σi表示矩阵S中保留的奇异值；
计算矩阵S中奇异值的均值，将该均值作为截断奇异值的阈值，其中，矩阵S中奇异值大于截断奇异值的阈值的保留，小于截断奇异值的阈值的进行零化处理；
则公式(7)的截断奇异值解
Figure BDA0002979540070000092
为：
Figure BDA0002979540070000093
公式(11)中，
Figure BDA0002979540070000094
根据最小二乘原理VTPV＝min，公式(8)写成以下形式：
(BTPB)X＝(BTPL) (12)
将公式(12)的等式左端和右端同时加未知参数KX，并化简得到公式(13)：
(BTPB+KI)X＝BTPL+KX (13)
由公式(13)整理得，最小二乘谱修正迭代公式为：
X＝(BTPB+KI)-1(BTPL+KX) (14)
公式(14)中，K为任意实数；
基于最小二乘谱修正迭代计算未知参数X的解，转到步骤II.3.5；
II.3.4.利用最小二乘法求解未知参数X的解，转到步骤II.3.5；
II.3.5.根据计算得到的未知参数X的解，得到未知参数X中包含的参数值，即测站的位置参数(dX,dY,dZ)、接收机钟差dtr、模糊度参数dN以及球谐展开的系数Cnm和Snm；
其中，将未知参数X中包含的位置参数(dX,dY,dZ)，分别加到测站的近似坐标(X0,Y0,Z0)上即精密单点定位求出的地心地固坐标系下的测站坐标；
II.3.6.将测站坐标由地心地固坐标系转换到站心坐标系，实现GNSS精密单点定位。
本发明具有如下优点：
如上所述，本发明述及了一种基于球谐展开的GNSS单点定位方法，该方法采用球谐展开描述与测站和卫星之间的高度角和方位角有关的误差项，包括对流层延迟误差以及电离层延迟误差，使用GNSS观测数据和卫星星历数据即可实现测站空间位置的计算。本发明能高效、快速地获得测点的位置信息，具有实际操作简单、数据处理便捷以及计算效率高等优势。
附图说明
图1为本发明实施例1中基于球谐展开的GNSS标准单点定位的流程示意图；
图2为本发明实施例2中基于球谐展开的GNSS精密单点定位的流程示意图；
图3为本发明实施例基于球谐展开的GNSS标准单点定位方法的结果统计图；
图4为本发明实施例基于球谐展开的GNSS精密单点定位方法的结果统计图。
具体实施方式
下面结合附图以及具体实施方式对本发明作进一步详细说明：
实施例1
本实施例述及了一种基于球谐展开的GNSS标准单点定位方法。如图1所示，该基于球谐展开的GNSS标准单点定位方法包括如下步骤：
I.1.建立利用伪距观测值进行的标准单点定位观测方程。
实现标准单点定位，读取观测数据和广播星历数据，处理伪距观测数据，计算与测站和卫星间的高度角和方位角无关或关系不密切的误差，包括卫星星历误差和卫星钟差。由于卫星星历误差的量级与标准单点定位误差的量级基本相同，因此，广播星历可用于低精度的标准单点定位，在标准单点定位中卫星星历误差可以暂时忽略。
在标准单点定位中，由广播星历内插计算得到每个历元对应的卫星位置；而卫星钟差使用广播星历提供的卫星钟差参数求得，一般卫星钟差ts由一个多项式拟合得到：ts＝α0+α1(ti-t0)+α2(ti-t0)2，α0、α1和α2分别为钟差、钟漂和钟的频漂系数，可以在导航电文中读取；t0为多项式计算卫星钟差的参考时刻，ti表示信号接收时刻即数据记录时刻。
地面接收机钟大多采用石英钟，地面接收机的精度低于卫星钟，钟差变化无规律性，无法通过多项式较好拟合，因此，将接收机钟差tr作为待求参数进行估计。
建立利用伪距观测值进行的标准单点定位观测方程，如公式(1)所示：
Figure BDA0002979540070000111
式中，i表示第i观测历元，j表示卫星编号；ρ表示伪距观测值，
Figure BDA0002979540070000112
表示卫星j在第i观测历元的伪距观测值；
Figure BDA0002979540070000113
表示接收机的位置与卫星j在第i观测历元的位置之间的几何距离；c表示光速；tr表示测站接收机钟差，(tr)i表示第i观测历元测站的接收机钟差；ts表示卫星钟差，
Figure BDA0002979540070000114
表示卫星j在第i观测历元的钟差；
Figure BDA0002979540070000115
表示卫星j在第i观测历元的信号传播路径上的对流层延迟误差；
Figure BDA0002979540070000116
表示卫星j在第i观测历元的信号传播路径上的电离层延迟误差；ερ表示伪距观测数据残差；定义卫星j在第i观测历元观测瞬间的空间三维坐标为
Figure BDA0002979540070000117
测站的近似坐标为(X0,Y0,Z0)，则卫星到测站近似位置的几何距离
Figure BDA0002979540070000118
表示为：
Figure BDA0002979540070000119
GNSS伪距观测数据类型分为单频观测数据和双频观测数据。若进行标准单点定位的GNSS接收机为双频接收机，则电离层延迟误差通过双频消电离层组合来消除，组合形式如下：
Figure BDA00029795400700001110
ρ为伪距消电离层组合观测值；f1、f2为相应观测值的频率；ρ1、ρ2分别表示两种频率对应的伪距观测值。若观测数据为双频观测值时，则公式(1)中，电离层延迟误差的一阶主项通过双频消电离层组合消除，电离层延迟误差项
Figure BDA00029795400700001111
为0。此时，本实施例1中的标准单点定位观测方程中，球谐展开主要用于表示对流层延迟误差项，球谐展开的阶次与单频伪距观测值构建的观测方程有所不同。
I.2.在标准单点定位中，对流层延迟误差和电离层延迟误差都是卫星信号穿过大气层时产生的一种信号延迟误差，都与卫星信号的传播有关；而卫星信号的传播路径又与测站和卫星的位置有关，利用测站坐标和卫星位置可以计算出测站和卫星之间的高度角和方位角，应用球谐展开，就能够表示与测站和卫星间的高度角和方位角有关的误差项，包括对流层延迟误差和电离层延迟误差，以消除或最大程度的削弱它们对标准单点定位结果的影响。基于球谐展开的标准单点定位观测方程表示为：
Figure BDA0002979540070000121
式中，n为球谐展开的阶数，m为球谐展开的次数，Nmax为球谐展开的最大阶数；
Figure BDA0002979540070000122
表示n阶m次的缔合勒让德多项式；Cnm和Snm分别表示n阶m次的球谐函数模型的系数，Cnm和Snm为球谐展开的待求参数；
Figure BDA0002979540070000123
表示测站与卫星j在第i观测历元之间的高度角，
Figure BDA0002979540070000124
表示测站与卫星j在第i观测历元之间的方位角。
经过公式(2)表示的基于球谐展开的标准单点定位观测方程，不涉及对流层延迟误差的改正模型，无需考虑不同模型在不同地区的适用问题，并且本实施例将球谐展开的系数当作待求参数直接求解，来改正与测站和卫星之间的高度角和方位角有关的误差项，主要包括对流层延迟误差，无需考虑气象信息，更加具有普适性。为了方便表示球谐函数展开，令：
Figure BDA0002979540070000125
则公式(2)简化为：
Figure BDA0002979540070000126
由简化后的标准单点定位观测方程(3)，得到标准单点定位误差方程(4)，公式如下：
Figure BDA0002979540070000127
其中，vρ表示伪距观测数据的改正值；
Figure BDA0002979540070000128
表示卫星j在第i观测历元的伪距观测数据的改正值。
I.3.基于公式(4)得到标准单点定位误差方程的线性化表达式，并进行简化处理，对简化后的标准单点定位误差方程的线性化表达式进行滑动解算。
I.3.1.对标准单点定位误差方程(4)在测站的近似坐标(X0,Y0,Z0)处进行泰勒级数展开并保留一阶项，则标准单点定位误差方程的线性化表达式为：
Figure BDA0002979540070000129
公式(5)中，令
Figure BDA00029795400700001210
Figure BDA0002979540070000131
Figure BDA0002979540070000132
Figure BDA0002979540070000133
表示测站近似坐标与卫星j在第i观测历元坐标计算出的待求参数dX的系数，待求参数dX是测站近似坐标X0的改正数；
Figure BDA0002979540070000134
表示测站近似坐标与卫星j在第i观测历元坐标计算出的待求参数dY的系数，待求参数dY是测站近似坐标Y0的改正数；
Figure BDA0002979540070000135
表示测站近似坐标与卫星j在第i观测历元坐标计算出的待求参数dZ的系数，待求参数dZ是测站近似坐标Z0的改正数；
Figure BDA0002979540070000136
表示由测站近似坐标与卫星j在第i观测历元坐标计算出的球谐展开的待求参数Cnm的系数，
Figure BDA0002979540070000137
表示由测站近似坐标与卫星j在第i观测历元坐标计算出的球谐展开的待求参数Snm的系数；tr表示测站接收机钟差，因为测站接收机钟差无法很好地被多项式拟合，所以将其作为待求参数进行解算，作为待求参数求解时，dtr的系数为1；d(tr)i表示在第i观测历元测站接收机钟差的待求参数。
简化后的标准单点定位误差方程的线性化表达式，如公式(6)所示；
Figure BDA0002979540070000138
公式(6)中，i＝1,2,…,epoch，epoch表示最大观测历元数；
Figure BDA0002979540070000139
表示由卫星j在第i观测历元伪距观测值、卫星j在第i观测历元的钟差以及测站近似坐标与卫星j在第i观测历元之间的几何距离计算出的常数项，
Figure BDA00029795400700001310
由公式(6)能够看出，此时，有epoch×j个观测方程。而待求参数包括3个位置参数(dX,dY,dZ)，epoch个接收机钟差dtr参数以及(Nmax+1)2个球谐展开的参数Cnm、Snm。
当epoch充分多时，能够通过最小二乘求得所有未知数的最优解。因此只要改正项是可以线性描述的，就能通过增加观测历元的方式来增加观测方程数量，从而解算出这些未知参数。
通过以上分析，为了确保未知参数的估计为最优解，本实施例选择合适的滑动窗口，对简化后的标准单点定位误差方程的线性表达式(6)进行滑动解算。
设定滑动窗口选择i个历元，每个历元最多能够观测到j颗卫星，则在这个滑动窗口下所有可视卫星的伪距观测值组成imax个方程，其中，imax＝i×j。
定义标准单点定位误差方程的系数矩阵、伪距观测值的改正数向量矩阵、常数项矩阵以及未知参数矩阵分别为矩阵B、V、L、X，则它们的表达式分别如下：
Figure BDA0002979540070000141
Figure BDA0002979540070000142
Figure BDA0002979540070000143
X＝[dX dY dZ (dtr)1 … (dtr)i C00 … CNmaxNmax SNmaxNmax]T。
则标准单点定位误差方程的线性化表达式(6)用一个矩阵形式表达，如公式(7)所示：
V＝BX-L (7)
在构建公式(7)的系数矩阵B时，只需要通过简单的数学计算，计算每一个待求参数(位置参数、测站接收机钟差参数及球谐展开的参数)的系数即可，计算过程简单，构建误差方程的速度迅速，实现标准单点定位的效率高。目前，与测站和卫星之间的高度角和方位角有关的误差项，如对流层延迟误差等还无法用模型精确改正，参考处理测站接收机钟差参数的方式，将其当做待求参数进行求解，是一种很好地方式。基于球谐展开的GNSS标准单点定位方法的提出正是基于这一思想，利用球谐展开表示与测站和卫星之间的高度角和方位角有关的误差项，在标准单点定位方程中直接求解球谐展开的系数。
未知参数X的最小二乘估值为：X＝(BTPB)-1BTPL (8)
公式(8)中，P为单位权矩阵。
I.3.2.如果标准单点定位误差方程的系数矩阵B是良态，则用公式(8)求出的未知参数的解是最优无偏估计，当系数矩阵B是病态时，则利用最小二乘法求出的不再是最优解。
因此，需要首先判断标准单点定位误差方程的系数矩阵B是否为病态，若系数矩阵B是病态，则执行步骤I.3.3；否则，若系数矩阵B为良态，执行步骤I.3.4。
本实施例通过计算的系数阵B的条件数，设置条件数的经验阈值，并将条件数与经验阈值作比较，当条件数大于经验阈值时，则表明误差方程(7)为病态，否则为良态。
I.3.3.为了解决误差方程的病态问题，本实施例提出了利用截断奇异值分解的方法为最小二乘谱修正迭代提供待求参数X的初值，基于最小二乘谱修正迭代计算待求参数X的值。截断奇异值法的思想是计算矩阵S中奇异值的均值，将该均值作为截断奇异值的阈值；筛选出小于阈值的奇异值将其零化。由于截断奇异值分解方法是将方程的最小二乘法转换为直接解法，因而避免了病态方程的解过分偏离真实解，有效解决了误差方程系数矩阵B的病态问题。
设系数矩阵B∈Rn×m，Rn×m表示n行m列的实数阵；则系数矩阵B的奇异值分解为：
B＝USVT (9)
公式(9)中，U∈n×n，V∈m×m，U、V均为正交矩阵，S∈n×m为一对角矩阵。
系数矩阵B∈Rn×m的截断奇异值矩阵Bk定义为：
Figure BDA0002979540070000151
公式(10)中，矩阵S中最小的(r-k)个非奇异值用零代替，即被截断，且k≤r；r表示系数矩阵B的秩，k表示矩阵S中保留的奇异值数量；ui表示矩阵U对应的向量，vi表示矩阵V对应的向量；σi表示矩阵S中保留的奇异值。计算矩阵S中奇异值的均值，将该值作为截断奇异值的阈值；矩阵S中奇异值大于该阈值的保留，小于该阈值的进行零化处理，则公式(7)的截断奇异值解
Figure BDA0002979540070000152
为：
Figure BDA0002979540070000153
公式(11)中，
Figure BDA0002979540070000154
根据最小二乘原理VTPV＝min，公式(8)写成以下形式：
(BTPB)X＝(BTPL) (12)
将公式(12)的等式左端和右端同时加未知参数KX，并化简得到公式(13)：
(BTPB+KI)X＝BTPL+KX (13)
由公式(13)整理得，最小二乘谱修正迭代公式为：
X＝(BTPB+KI)-1(BTPL+KX) (14)
公式(14)中，K为任意实数。最小二乘谱修正迭代是一种基于最小二乘的迭代计算方法；基于最小二乘谱修正迭代计算未知参数X的解。求解过程如下：
基于最小二乘谱修正迭代解算X包括迭代①和迭代②两个迭代过程；其中，设定用于判断迭代①是否收敛的第一比较阈值，设定用于判断迭代②是否收敛的第二比较阈值。迭代①第一次迭代时，设定K的初始值为1，将公式(11)得到的截断奇异值解
Figure BDA0002979540070000161
作为公式(14)第一次迭代时X的初值，并赋值给公式(14)等式右边的X；
每次迭代过程中，将上一次迭代时利用公式(14)得到的未知参数X的值，作为本次迭代过程中X的初值，并赋值给公式(14)等式右边的X；
基于公式(14)求解得到每次迭代过程中未知参数X的值；
将每次迭代过程中求解得到的X的值，与每次迭代过程中X的初值进行差值运算；
若经过差值运算，得到的差值大于第一比较阈值，则表示迭代不收敛，需要改正最小二乘谱修正迭代中的K值，即将K值递增1，继续上述迭代①过程；
若经过差值运算得到的差值小于或等于第一比较阈值，则跳出迭代①，进入迭代②过程。
迭代②将利用公式(14)求解得到的未知参数X的值与第二比较阈值进行比较；
若利用公式(14)求解得到的未知参数X的值大于第二比较阈值，则表示迭代②未收敛，此时，将位置参数(dX,dY,dZ)分别加到测站的近似坐标(X0,Y0,Z0)上，更新测站的近似坐标，通过公式(5)对更新测站近似坐标后的标准单点定位误差方程进行线性化，通过公式(6)简化线性化后的标准单点定位误差方程；重新进入迭代①过程求解出未知参数X的值。
若利用公式(14)求解得到的未知参数X的值小于或等于第二比较阈值，则表示迭代②收敛，跳出迭代；此时利用公式(14)得到的未知参数X的值，即未知参数X的解。
求得未知参数X的解后，转到步骤I.3.5。
I.3.4.利用最小二乘法求解未知参数X的解，转到步骤I.3.5。
步骤I.3.4中，基于最小二乘法求解未知参数X的解的过程如下：
设定用于判断迭代是否收敛的第三比较阈值；每次迭代过程中，利用公式(8)求解出未知参数X的值，将本次迭代得到的未知参数X的值与第三比较阈值进行比较；
若利用公式(8)求解得到的未知参数X的值大于第三比较阈值，则表示迭代未收敛，此时，将位置参数(dX,dY,dZ)，分别加到测站的近似坐标(X0,Y0,Z0)上，更新测站的近似坐标，执行公式(5)-公式(8)的计算过程，并重新求解未知参数X的值；
若利用公式(8)求解得到的未知参数X的值小于或等于第三比较阈值，则表示迭代收敛，跳出迭代；此时，利用公式(8)求解得到的未知参数X的值，即未知参数X的解。
I.3.5.根据计算得到的未知参数X的解，得到未知参数X中包含的参数值，即测站的位置参数(dX,dY,dZ)、接收机钟差dtr以及球谐展开的系数Cnm和Snm；
其中，未知参数X中的位置参数(dX,dY,dZ)分别加到测站的近似坐标(X0,Y0,Z0)上，即标准单点定位求出的地心地固坐标系下的测站坐标(X,Y,Z)。
通过求解得到的球谐展开的系数Cnm和Snm可直接用于后续标准单点定位。
具体应用过程为：读取已有的球谐展开系数，通过球谐展开表示与测站和卫星之间的高度角和方位角有关的误差项，主要包括对流层延迟误差和电离层延迟误差；构建单点定位观测方程时，不再考虑与测站和卫星之间的高度角和方位角有关的误差项的改正问题，省略定位中对流层延迟误差和电离层延迟误差的求解过程，使得定位方程仅仅需要解算测站的位置参数、测站接收机钟差参数即可，极大地简化了标准单点定位流程。
I.3.6.将测站坐标由地心地固坐标系转换到站心坐标系，实现GNSS标准单点定位。
将测站坐标由地心地固坐标系转换到站心坐标系的过程如下：
地心地固坐标系转经纬高坐标系的公式为：
Figure BDA0002979540070000171
式中，(X,Y,Z)为在地心地固坐标系的坐标，(lon,lat,alt)为在经纬高坐标系的坐标；M表示基准椭球体的曲率半径，E表示椭球偏心率；根据测站近似坐标(X0,Y0,Z0)，利用公式(15)计算得到该测站近似坐标(X0,Y0,Z0)在经纬高坐标系下的对应坐标，即(lon0,lat0,alt0)。
地心地固坐标系转换到站心坐标系的公式为：
Figure BDA0002979540070000172
式中，(e1,n1,u1)为以测站近似坐标(X0,Y0,Z0)为原点求出的(X,Y,Z)坐标的位置。
S为坐标转换矩阵，
Figure BDA0002979540070000173
根据以上坐标转换关系，将测站坐标由地心地固坐标系转换到站心坐标系。
本实施例1中的GNSS标准单点定位方法，基于球谐展开表示与测站和卫星间的方位角和高度角有关的误差项，从而能够快速、高效、准确地进行标准单点定位。
实施例2
本实施例2述及了一种基于球谐展开的GNSS精密单点定位方法。如图2所示，该基于球谐展开的GNSS精密单点定位方法包括如下步骤：
II.1.建立利用载波相位观测值进行的精密单点定位观测方程。
读取观测数据、广播星历数据以及精密星历数据，对载波相位观测数据进行预处理，计算与测站和卫星间的高度角和方位角无关或关系不密切的误差；以上误差主要包括计算卫星钟差、探测周跳以及计算载波相位整周模糊度参数的初始值等。
对于精密单点定位，广播星历的精度无法满足精密单点定位的要求，因此必须采用卫星精密星历计算卫星位置；利用二阶多项式估计的卫星钟差无法满足精密单点定位的精度要求；由于不同的卫星钟差均不相同，必须事先确定其大小，再代入到观测方程消除其影响。
目前，IGS提供的15min、5min采样间隔的事后精密星历和5min、30s采样间隔的卫星钟差产品，能够满足精密定位的精度要求。观测值的采样间隔一般都小于上述文件中数据的采样间隔，因此，需要内插计算得到每个历元对应的卫星位置和卫星钟差。
地面接收机钟大多采用石英钟，精度低于卫星钟，钟差变化无规律性，无法通过多项式较好拟合，因此，将接收机钟差tr作为待求参数进行估计。
建立的利用载波相位观测值进行的精密单点定位观测方程，如公式(1)所示。
Figure BDA0002979540070000181
公式(1)中，i表示第i观测历元；j表示卫星编号；
Figure BDA0002979540070000182
表示载波相位观测值，将载波相位原始观测值乘上相应的波长λ，得到以距离表示的载波相位观测数据，称为等效距离；
Figure BDA0002979540070000183
表示卫星j在第i观测历元的载波相位观测值，即等效距离；
Figure BDA0002979540070000184
表示接收机的位置与卫星j在第i观测历元位置之间的几何距离；tr表示接收机钟差，(tr)i表示第i观测历元接收机钟差；ts表示卫星钟差，
Figure BDA0002979540070000185
表示卫星j在第i观测历元的钟差；
Figure BDA0002979540070000186
表示卫星j在第i观测历元的信号传播路径上的对流层延迟误差；
Figure BDA0002979540070000187
表示卫星j在第i观测历元的信号传播路径上的电离层延迟误差；c表示光速；λ表示载波的波长；N表示卫星相位观测数据的整周模糊度参数，
Figure BDA0002979540070000188
表示卫星j在第i观测历元的相位观测数据整周模糊度参数；
Figure BDA0002979540070000189
表示载波相位观测数据残差。定义卫星j在第i历元观测瞬间的空间三维坐标为
Figure BDA00029795400700001810
测站的近似坐标为(X0,Y0,Z0)，则卫星到测站近似位置的几何距离
Figure BDA00029795400700001811
表示为：
Figure BDA00029795400700001812
GNSS载波相位观测数据类型分为单频观测数据和双频观测数据。
若进行标准单点定位的GNSS接收机为双频接收机，则电离层延迟误差可以通过双频消电离层组合来消除，具体的组合形式如下：
Figure BDA0002979540070000191
式中，
Figure BDA0002979540070000192
表示载波相位消电离层组合观测值(等效距离)；f1、f2表示相应观测值的频率；
Figure BDA0002979540070000193
分别表示两种频率对应的载波相位观测值(等效距离)。若观测数据为双频观测值时，则上述公式(1)中，电离层延迟误差的一阶主项通过双频消电离层组合消除，电离层延迟误差项
Figure BDA0002979540070000194
为0。
此时，本实施例2中的精密单点定位观测方程中，球谐展开主要用于表示对流层延迟误差项，球谐展开的阶次与单频载波相位观测值构建的观测方程有所不同。
在精密单点定位中，利用载波相位观测值组成GF(几何无关)组合观测值，进行整周模糊度的解算和周跳的探测。具体的组合形式如下：
Figure BDA0002979540070000195
式中，
Figure BDA0002979540070000196
表示GF组合载波相位观测值(等效距离)；
Figure BDA0002979540070000197
分别表示两种频率对应的载波相位观测值；λ1、λ2分别表示两种频率对应的载波相位的波长；f1、f2表示两种载波相位的频率。GF组合观测值与几何距离无关，因此，该组合适合用来探测周跳。
II.2.基于球谐展开表示与测站和卫星间的高度角和方位角有关的误差项，其中，误差项主要包括对流层延迟误差以及电离层延迟误差。在精密单点定位中，对流层延迟误差和电离层延迟误差都是卫星信号穿过大气层时产生的一种信号延迟误差，都与卫星信号的传播有关；而卫星信号的传播路径又与测站和卫星的位置有关，利用测站坐标和卫星位置可以计算出测站和卫星之间的高度角和方位角，应用球谐展开，就可以表示与测站和卫星间的高度角和方位角有关的误差项，主要包括对流层延迟误差和电离层延迟误差，以消除或最大程度的削弱其对精密单点定位结果的影响。经过球谐展开后的精密单点定位观测方程表示为：
Figure BDA0002979540070000198
式(2)中，n为球谐展开的阶数，m为球谐展开的次数，Nmax为球谐展开的最大阶数；Pnm(cose)表示n阶m次的缔合勒让德多项式；Cnm和Snm表示n阶m次的球谐展开的系数；
Figure BDA0002979540070000199
表示测站与卫星j在第i观测历元之间的高度角，
Figure BDA00029795400700001910
表示测站与卫星j在第i观测历元之间的方位角。公式(2)提出的基于球谐展开的精密单点定位观测方程，不涉及对流层延迟误差的改正模型，无需考虑不同模型在不同地区的适用问题，并且本实施例将球谐展开的系数当作待求参数直接求解，来改正与测站和卫星之间的高度角和方位角有关的误差项，主要包括对流层延迟误差，无需考虑气象信息，更加具有普适性；为了方便表示球谐函数展开，令：
Figure BDA0002979540070000201
则公式(2)简化为：
Figure BDA0002979540070000202
由公式(3)中的精密单点定位观测方程，得到精密单点定位误差方程：
Figure BDA0002979540070000203
其中，
Figure BDA0002979540070000204
表示载波相位观测数据(等效距离)的改正值。
Figure BDA0002979540070000205
表示第i观测历元卫星j的载波相位观测数据(等效距离)的改正值。
II.3.基于公式(4)得到精密单点定位误差方程的线性化表达式，并简化处理，然后对简化后的精密单点定位误差方程的线性化表达式进行滑动解算。
II.3.1.对精密单点定位误差方程(4)在测站的近似坐标(X0,Y0,Z0)处进行泰勒级数展开并保留一阶项，则精密单点定位误差方程的线性化表达式为：
Figure BDA0002979540070000206
公式(5)中，
Figure BDA0002979540070000207
表示卫星j在第i历元的模糊度初值；
Figure BDA0002979540070000208
表示卫星j在第i历元的相位观测数据整周模糊度参数，dN表示相位观测数据整周模糊度参数的改正数，
Figure BDA0002979540070000209
表示卫星j在第i历元的相位观测数据整周模糊度参数的改正数，在精密单点定位中，由于整周模糊度参数不能被很好地模型化，因此将其改正数作为待求参数进行求解，作为待求参数求解时，
Figure BDA00029795400700002011
的系数为1；d(tr)i表示在第i观测历元测站接收机钟差的待求参数，与整周模糊度的处理方式相似，因为测站接收机钟差无法很好地被多项式拟合，因此将其作为待求参数进行求解，作为待求参数求解时，dtr的系数为1。
令
Figure BDA00029795400700002010
Figure BDA0002979540070000211
Figure BDA0002979540070000212
Figure BDA0002979540070000213
为测站近似坐标与卫星j在第i观测历元坐标计算出的待求参数dX的系数，待求参数dX是测站近似坐标X0的改正数；
Figure BDA0002979540070000214
为测站近似坐标与卫星j在第i观测历元坐标计算出的待求参数dY的系数，待求参数dY是测站近似坐标Y0的改正数；
Figure BDA0002979540070000215
为测站近似坐标与卫星j在第i观测历元坐标计算出的待求参数dZ的系数，待求参数dZ是测站近似坐标Z0的改正数；Anm表示球谐展开待求参数Cnm的系数，Bnm表示球谐展开待求参数Snm的系数；
Figure BDA0002979540070000216
表示由测站近似坐标与卫星j在第i观测历元坐标计算出的球谐展开的待求参数Cnm的系数，
Figure BDA0002979540070000217
为由测站近似坐标与卫星j在第i观测历元坐标计算出的球谐展开的待求参数Snm的系数。
简化后的精密单点定位误差方程的线性化表达式，如公式(6)所示：
Figure BDA0002979540070000218
公式(6)中，i＝1,2,,…,epoch，epoch表示最大观测历元；
Figure BDA0002979540070000219
表示由卫星j在第i观测历元载波相位观测值、卫星j在第i观测历元的钟差、卫星j在第i历元的相位观测数据整周模糊度初始值以及测站近似坐标与卫星j在第i观测历元之间的几何距离计算出的常数项，
Figure BDA00029795400700002110
由公式(6)可以看出，此时有epoch×j个观测方程。而待求参数包括：3个位置参数(dX,dY,dZ)，epoch个接收机钟差参数dtr，j个模糊度参数dN以及(Nmax+1)2个球谐展开的系数Cnm、Snm。
当epoch充分多时，能够通过最小二乘求得所有未知数的最优解，只要改正项是可以线性描述的，就可以通过增加观测历元的方式来增加观测方程数量，从而解算出这些未知参数。
设定滑动窗口选择i个历元，每个历元最多能够观测到j颗卫星，则在这个滑动窗口下所有可视卫星的载波相位观测值组成imax个方程，其中，imax＝i×j。
定义精密单点定位误差方程的系数矩阵、载波相位观测值的改正数向量矩阵、常数项矩阵以及未知参数矩阵分别为矩阵B、V、L、X，则它们的表达式分别如下：
Figure BDA0002979540070000221
Figure BDA0002979540070000222
Figure BDA0002979540070000223
Figure BDA0002979540070000224
Figure BDA0002979540070000225
表示第1颗卫星在第1个观测历元的整周模糊度参数，
Figure BDA0002979540070000226
表示第j颗卫星在第i个观测历元的整周模糊度参数，精密单点定位中求得的是整周模糊度初始值的改正数。
则精密单点定位误差方程的线性化表达式(6)用一个矩阵形式表达，如公式(7)所示；
V＝BX-L (7)
构建精密单点定位误差方程式(7)的系数矩阵B时，只需要通过简单的数学计算，计算每一个待求参数(位置参数、测站接收机钟差参数、整周模糊度参数及球谐展开的参数)的系数即可，计算过程简单，构建误差方程的速度迅速，实现精密单点定位的效率高。目前，与测站和卫星之间的高度角和方位角有关的误差项，如对流层延迟误差等还无法用模型精确改正，参考处理测站接收机钟差参数的方式，将其当做待求参数进行求解，是一种很好地方式；基于球谐展开的GNSS精密单点定位方法的提出正是基于这一思想，利用球谐展开表示与测站和卫星之间的高度角和方位角有关的误差项，在精密单点定位方程中直接求解球谐展开的系数。未知参数X的最小二乘估值为：X＝(BTPB)-1BTPL (8)
公式(8)中，P为单位权矩阵。
II.3.2.如果精密单点定位误差方程的系数矩阵B是良态，则用公式(8)求出的未知参数的解时最优无偏估计，当系数矩阵B是病态时，则利用最小二乘法求出的不再是最优解。
因此，需要首先判断精密单点定位误差方程的系数矩阵B是否为病态，若系数矩阵B是病态，则执行步骤II.3.3；否则，系数矩阵B为良态，执行步骤II.3.4。
本实施例通过计算精密单点定位误差方程的系数阵B的条件数，设置条件数的经验阈值，将条件数与经验阈值作比较，当条件数大于经验阈值时，则误差方程为病态。
II.3.3.为了解决误差方程的病态问题，本实施例提出了利用截断奇异值分解的方法，为最小二乘谱修正迭代提供待求参数X的初值，基于最小二乘谱修正迭代计算待求参数X的值。由于截断奇异值分解方法是将方程的最小二乘法转换为直接解法，因而避免了病态方程的解过分偏离真实解，有效解决了误差方程系数矩阵B的病态问题。
设系数矩阵B∈Rn×m，Rn×m表示n行m列的实数阵，则系数矩阵B的奇异值分解为：
B＝USVT (9)
公式(9)中，U∈n×n，V∈m×m，U、V均为正交矩阵，S∈n×m为一对角矩阵。
系数矩阵B∈Rn×m的截断奇异值矩阵Bk定义为：
Figure BDA0002979540070000231
公式(10)中，矩阵S中最小的(r-k)个非奇异值用零代替，即被截断，且k≤r；r表示系数矩阵B的秩，k表示矩阵U中保留的奇异值数量，ui表示矩阵U对应的向量，vi表示矩阵V对应的向量；σi表示矩阵U中保留的奇异值。计算矩阵S中奇异值的均值，将该值作为截断奇异值的阈值；矩阵S中奇异值大于该截断奇异值的阈值的保留，小于该截断奇异值的阈值的进行零化处理。则公式(7)的截断奇异值解
Figure BDA0002979540070000232
为：
Figure BDA0002979540070000233
公式(11)中，
Figure BDA0002979540070000234
最小二乘谱修正迭代是一种基于最小二乘的迭代计算方法。基于最小二乘谱修正迭代计算待求参数X；根据最小二乘原理VTPV＝min，公式(8)写成以下形式：
(BTPB)X＝(BTPL) (12)
将公式(12)的等式左端和右端同时加未知参数KX，并化简得到公式(13)：
(BTPB+KI)X＝BTPL+KX (13)
由公式(13)整理得，最小二乘谱修正迭代公式为：
X＝(BTPB+KI)-1(BTPL+KX) (14)
公式(14)中，K为任意实数；基于最小二乘谱修正迭代计算未知参数X的解。
基于最小二乘谱修正迭代公式计算X的解的过程如下：
基于最小二乘谱修正迭代解算X包括迭代①和迭代②两个迭代过程；其中，设定用于判断迭代①是否收敛的第一比较阈值，设定用于判断迭代②是否收敛的第二比较阈值；
迭代①第一次迭代时，设定K的初始值为1，将公式(11)得到的截断奇异值解
Figure BDA0002979540070000241
作为公式(14)第一次迭代时X的初值，并赋值给公式(14)等式右边的X；
每次迭代过程中，将上一次迭代时利用公式(14)得到的未知参数X的值，作为本次迭代过程中X的初值，并赋值给公式(14)等式右边的X；
基于公式(14)求解得到每次迭代过程中未知参数X的值；
将每次迭代过程中求解得到的X的值，与每次迭代过程中X的初值进行差值运算；
若经过差值运算，得到的差值大于第一比较阈值，则表示迭代不收敛，需要改正最小二乘谱修正迭代中的K值，即将K值递增1，继续上述迭代①过程；
若经过差值运算得到的差值小于或等于第一比较阈值，则跳出迭代①，进入迭代②过程；
迭代②将利用公式(14)求解得到的未知参数X的值与第二比较阈值进行比较；
若利用公式(14)求解得到的未知参数X的值大于第二比较阈值，则表示迭代②未收敛，此时，将位置参数(dX,dY,dZ)分别加到测站的近似坐标(X0,Y0,Z0)上，更新测站近似坐标，通过公式(5)对更新测站近似坐标后的精密单点定位误差方程进行线性化，通过公式(6)简化线性化后的精密单点定位误差方程；重新进入迭代①过程求解出未知参数X的值；
若利用公式(14)求解得到的未知参数X的值小于或等于第二比较阈值，则表示迭代②收敛，跳出迭代；此时利用公式(14)得到的未知参数X的值，即待求解的未知参数X的解。
未知参数X的解后，转到步骤II.3.5。
II.3.4.利用最小二乘法求解未知参数X的解，转到步骤II.3.5。
步骤II.3.4中，基于最小二乘法求解未知参数X的解的过程如下：
设定用于判断迭代是否收敛的第三比较阈值；每次迭代过程中，利用公式(8)求解出未知参数X的值，将本次迭代得到的未知参数X的值与第三比较阈值进行比较；
若利用公式(8)求解得到的未知参数X的值大于第三比较阈值，则表示迭代未收敛，此时，将位置参数(dX,dY,dZ)分别加到测站的近似坐标(X0,Y0,Z0)上，更新测站近似坐标；
执行公式(5)-公式(8)的计算过程，并重新求解未知参数X的值；
若利用公式(8)求解得到的未知参数X的值小于或等于第三比较阈值，则表示迭代收敛，跳出迭代；此时，利用公式(8)求解得到的未知参数X的值，即待求解的未知参数X的解。
II.3.5.根据计算得到的未知参数X的解，得到未知参数X中包含的参数值，即测站的位置参数(dX,dY,dZ)、接收机钟差dtr、模糊度参数dN以及球谐展开的系数Cnm和Snm；
其中，位置参数(dX,dY,dZ)分别加到测站的近似坐标(X0,Y0,Z0)上，即精密单点定位求出的地心地固坐标系下的测站坐标(X,Y,Z)。通过求解得到的球谐展开的系数Cnm和Snm，可以直接用于后续精密单点定位，具体应用过程如下：
读取已有的球谐展开系数，通过球谐展开表示与测站和卫星之间的高度角和方位角有关的误差项，主要包括对流层延迟误差和电离层延迟误差，构建的精密单点定位观测方程可以不再考虑与测站和卫星之间的高度角和方位角有关的误差项的改正问题，省略定位中的对流层延迟误差和电离层延迟误差的求解过程，这将使得定位方程仅仅需要解算测站的位置参数、测站接收机钟差参数以及整周模糊度参数即可，极大地简化了精密单点定位流程。
II.3.6.将测站坐标由地心地固坐标系转换到站心坐标系，实现GNSS精密单点定位。其中，本实施例2中坐标转换过程与上述实施例1中坐标转换过程相同，此处不再赘述。
本实施例2中述及的GNSS精密单点定位方法，基于球谐展开表示与测站和卫星间的方位角和高度角有关的误差项，从而能够快速、高效、准确地进行精密单点定位。
下面对利用本发明中基于球谐展开的GNSS单点定位方法所能达到的精度进行说明：
以上海佘山站为例：
获取2019年1月7日佘山站的GNSS观测数据和卫星星历数据。
GNSS观测数据采样间隔为30s，观测数据选择GPS卫星的双频伪距观测值以及双频载波相位观测值，站点的参考坐标为IGS提供的站点坐标。
1.首先利用本发明中基于球谐展开的GNSS标准单点定位方法获取站点位置信息。
观测数据选择GPS卫星双频伪距观测值，则电离层延迟误差通过伪距线性组合消除。
通过球谐展开表示与卫星和测站间的高度角和方位角有关的误差项，主要是对流层延迟误差。球谐展开各阶次的系数、以及测站的接收机钟差需要当作待求参数求解。
将球谐函数展开到3阶，滑动窗口选择10个历元。待求参数包括：点位坐标3个方向的改正值，10个接收机钟差(每个历元一个接收机钟差参数)，球谐展开到3阶3次的16个系数：A00、A10、A11、B11、A20、A21、B21、A22、B22、A30、A31、B31、A32、B32、A33、B33。
卫星的截止高度角设置为5°，最小二乘谱修正迭代的系数k设为2。
GNSS标准单点定位方法的定位结果如图3所示，统计结果见表1。
表1标准单点定位精度(单位/m)
E N U
MAX 1.2662 1.7738 4.9711
MIN -1.3027 -2.0436 -5.0853
MEAN 0.1101 -0.0136 0.3269
STD 0.3727 0.3952 1.3887
RMS 0.3886 0.3954 1.4267
表2 Bernese定位精度(单位/m)
E N U
MAX 6.0917 5.7515 9.7484
MIN -3.8416 -5.8132 -11.1756
MEAN 0.1268 -0.1259 0.5520
STD 1.2744 1.6002 3.0886
RMS 1.2807 1.6052 3.1375
通过表1统计结果不难发现，本发明标准单点定位方法得到的定位结果在E方向上的RMS优于0.39m，在N方向上的RMS优于0.40m，在U方向上的RMS优于1.43m。
由瑞士伯尔尼大学天文研究所开发的国际知名GNSS高精度数据处理软件Bernese软件解算的标准单点定位结果，如表2所示。
可见，GNSS高精度数据处理软件Bernese软件解算的定位结果在E方向上的RMS优于1.28m，在N方向上的RMS优于1.61m，在U方向上的RMS优于3.14m。
从统计表中看出，本发明标准单点定位结果与Bernese软件标准单点定位结果精度相当。
2.利用本发明中基于球谐展开的GNSS标准单点定位方法获取站点位置信息。
观测数据选择GPS卫星双频载波相位观测值，则电离层延迟误差通过载波相位线性组合消除。通过球谐展开表示与卫星和测站间的高度角和方位角有关的误差项，主要是对流层延迟误差；球谐展开各阶次的系数需要当作待求参数求解。
每颗卫星载波相位观测数据数据的整周模糊度以及测站的接收机钟差也需要作为待求参数求解。将球谐函数展开到6阶，滑动窗口选择15个历元。
待求参数包括：点位坐标3个方向的改正值，15个接收机钟差(每个历元一个接收机钟差参数)，每个历元下每颗卫星的模糊度参数，球谐展开到6阶6次的49个系数：A00、A10、A11、B11、A20、A21、B21、A22、B22、A30、A31、B31、A32、B32、A33、B33、A40、A41、B41、A42、B42、A43、B43、A44、B44、A50、A51、B51、A52、B52、A53、B53、A54、B54、A55、B55、A60、A61、B61、A62、B62、A63、B63、A64、B64、A65、B65、A66、B66。
卫星的截止高度角设置为5°，最小二乘谱修正迭代的系数k设为3。
GNSS精密单点定位方法的结果如图4所示，统计结果见表3。
表3精密单点定位精度(单位/m)
Figure BDA0002979540070000261
Figure BDA0002979540070000271
表4 Bernese定位精度(单位/m)
E N U
MAX 0.0274 0.0301 0.0461
MIN -0.0178 -0.0209 -0.0738
MEAN 0.0036 0.0048 -0.0082
STD 0.0075 0.0092 0.0200
RMS 0.0083 0.0104 0.0216
通过表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。
从统计表中看出，本发明精密单点定位结果与Bernese软件精密单点定位结果精度相当。
本发明基于球谐展开的GNSS单点定位方法，利用观测数据和卫星星历数据高效快速获得测点的点位信息，相比于现有的利用经验模型进行对流层延迟改正方法，具有解算效率高，无需提供气象文件，操作简便、工作难度小，测量结果可靠，测量精度高等优势。
当然，以上说明仅仅为本发明的较佳实施例，本发明并不限于列举上述实施例，应当说明的是，任何熟悉本领域的技术人员在本说明书的教导下，所做出的所有等同替代、明显变形形式，均落在本说明书的实质范围之内，理应受到本发明的保护。

Claims (4)

1.一种基于球谐展开的GNSS标准单点定位方法，其特征在于，包括如下步骤：
I.1.建立利用伪距观测值进行的标准单点定位观测方程，如公式(1)所示：
Figure FDA0002979540060000011
式中，i表示第i观测历元，j表示卫星编号；ρ表示伪距观测值，
Figure FDA0002979540060000012
表示卫星j在第i观测历元的伪距观测值；
Figure FDA0002979540060000013
表示接收机的位置与卫星j在第i观测历元的位置之间的几何距离；c表示光速；tr表示测站接收机钟差，(tr)i表示第i观测历元测站的接收机钟差；ts表示卫星钟差，
Figure FDA0002979540060000014
表示卫星j在第i观测历元的钟差；
Figure FDA0002979540060000015
表示卫星j在第i观测历元的信号传播路径上的对流层延迟误差；
Figure FDA0002979540060000016
表示卫星j在第i观测历元的信号传播路径上的电离层延迟误差；ερ表示伪距观测数据残差；
定义卫星j在第i观测历元观测瞬间的空间三维坐标为
Figure FDA0002979540060000017
测站的近似坐标为(X0,Y0,Z0)，则卫星到测站近似位置的几何距离
Figure FDA0002979540060000018
表示为：
Figure FDA0002979540060000019
I.2.基于球谐展开表示与测站和卫星之间的高度角和方位角有关的误差项，误差项包括对流层延迟误差以及电离层延迟误差；基于球谐展开的标准单点定位观测方程表示为：
Figure FDA00029795400600000110
式中，n为球谐展开的阶数，m为球谐展开的次数，Nmax为球谐展开的最大阶数；
Figure FDA00029795400600000111
表示n阶m次的缔合勒让德多项式；Cnm和Snm分别表示n阶m次的球谐展开的系数，Cnm和Snm为球谐展开的待求参数；
Figure FDA00029795400600000112
分别表示测站与卫星j在第i观测历元之间的高度角和方位角；为了方便表示球谐展开，令：
Figure FDA00029795400600000113
则公式(2)简化为：
Figure FDA00029795400600000114
由简化后的标准单点定位观测方程(3)，得到标准单点定位误差方程(4)，公式如下：
Figure FDA00029795400600000115
其中，vρ表示伪距观测数据的改正值；
Figure FDA0002979540060000021
表示卫星j在第i观测历元的伪距观测数据的改正值；
I.3.基于公式(4)得到标准单点定位误差方程的线性化表达式，并进行简化处理，然后对简化后的标准单点定位误差方程的线性化表达式进行滑动解算；
I.3.1.对标准单点定位误差方程(4)在测站的近似坐标(X0,Y0,Z0)处进行泰勒级数展开并保留一阶项，则标准单点定位误差方程的线性化表达式如公式(5)所示；
Figure FDA0002979540060000022
公式(5)中，令：
Figure FDA0002979540060000023
Figure FDA0002979540060000024
Figure FDA0002979540060000025
其中，
Figure FDA0002979540060000026
表示测站近似坐标与卫星j在第i观测历元坐标计算出的待求参数dX的系数，待求参数dX是测站近似坐标X0的改正数；
Figure FDA0002979540060000027
表示测站近似坐标与卫星j在第i观测历元坐标计算出的待求参数dY的系数，待求参数dY是测站近似坐标Y0的改正数；
Figure FDA0002979540060000028
表示测站近似坐标与卫星j在第i观测历元坐标计算出的待求参数dZ的系数，待求参数dZ是测站近似坐标Z0的改正数；
Figure FDA0002979540060000029
表示由测站近似坐标与卫星j在第i观测历元坐标计算出的球谐展开的待求参数Cnm的系数；
Figure FDA00029795400600000210
表示由测站近似坐标与卫星j在第i观测历元坐标计算出的球谐展开的待求参数Snm的系数；d(tr)i表示在第i观测历元测站接收机钟差的待求参数；
则简化后的标准单点定位误差方程的线性化表达式，如公式(6)所示；
Figure FDA00029795400600000211
公式(6)中，i＝1,2,…,epoch，epoch表示最大观测历元数；
Figure FDA00029795400600000212
表示由卫星j在第i观测历元伪距观测值、卫星j在第i观测历元的钟差以及测站近似坐标与卫星j在第i观测历元之间的几何距离计算出的常数项，其中，
Figure FDA0002979540060000031
设定滑动窗口选择i个历元，每个历元最多能够观测到j颗卫星，则在这个滑动窗口下所有可视卫星的伪距观测值组成imax个方程，其中，imax＝i×j；
定义标准单点定位误差方程的系数矩阵、伪距观测值的改正数向量矩阵、常数项矩阵以及未知参数矩阵分别为矩阵B、V、L、X，则它们的表达式分别如下：
Figure FDA0002979540060000032
Figure FDA0002979540060000033
Figure FDA0002979540060000034
X＝[dX dY dZ (dtr)1 … (dtr)i C00 … CNmaxNmax SNmaxNmax]T；
则标准单点定位误差方程的线性化表达式(6)用一个矩阵形式表达，如公式(7)所示：
V＝BX-L (7)
未知参数X的最小二乘估值为：X＝(BTPB)-1BTPL (8)
公式(8)中，P为单位权矩阵；
I.3.2.判断公式(7)中系数矩阵B是否为病态，经过判断，若系数矩阵B是病态，则执行步骤I.3.3；否则，系数矩阵B是良态，执行步骤I.3.4；
I.3.3.首先利用截断奇异值分解的方法，为最小二乘谱修正迭代提供待求解的未知参数X的初值，然后基于最小二乘谱修正迭代计算未知参数X的值；
设系数矩阵B∈Rn×m，Rn×m表示n行m列的实数阵；则系数矩阵B的奇异值分解为：
B＝USVT (9)
公式(9)中，U∈n×n，V∈m×m，U、V均为正交矩阵，S∈n×m为一对角矩阵；
系数矩阵B∈Rn×m的截断奇异值矩阵Bk定义为：
Figure FDA0002979540060000041
公式(10)中，矩阵S中最小的(r-k)个非奇异值用零代替，即被截断，且k≤r；r表示系数矩阵B的秩，k表示矩阵S中保留的奇异值数量；
ui表示矩阵U对应的向量，vi表示矩阵V对应的向量，σi表示矩阵S中保留的奇异值；
计算矩阵S中奇异值的均值，将该均值作为截断奇异值的阈值，其中，矩阵S中奇异值大于截断奇异值的阈值的保留，小于截断奇异值的阈值的进行零化处理；
则公式(7)的截断奇异值解
Figure FDA0002979540060000042
为：
Figure FDA0002979540060000043
公式(11)中，
Figure FDA0002979540060000044
根据最小二乘原理VTPV＝min，公式(8)写成以下形式：
(BTPB)X＝(BTPL) (12)
将公式(12)的等式左端和右端同时加未知参数KX，并化简得到公式(13)：
(BTPB+KI)X＝BTPL+KX (13)
由公式(13)整理得，最小二乘谱修正迭代公式为：
X＝(BTPB+KI)-1(BTPL+KX) (14)
公式(14)中，K为任意实数；
基于最小二乘谱修正迭代计算未知参数X的解，转到步骤I.3.5；
I.3.4.利用最小二乘法求解未知参数X的解，转到步骤I.3.5；
I.3.5.根据计算得到的未知参数X的解，得到未知参数X中包含的参数值，即测站的位置参数(dX,dY,dZ)、接收机钟差dtr以及球谐展开的系数Cnm和Snm；
将未知参数X中包含的位置参数(dX,dY,dZ)，分别加到测站的近似坐标(X0,Y0,Z0)上，即标准单点定位求出的地心地固坐标系下的测站坐标(X,Y,Z)；
I.3.6.将测站坐标由地心地固坐标系转换到站心坐标系，实现GNSS标准单点定位。
2.根据权利要求1所述基于球谐展开的GNSS标准单点定位方法，其特征在于，
所述步骤I.3.2中，系数矩阵B是否为病态的判断过程如下：计算公式(7)的系数阵B的条件数，设置条件数的经验阈值；将条件数与经验阈值进行比较，经过比较判断：
当条件数大于经验阈值时，则表明系数矩阵B为病态，否则，系数矩阵B为良态。
3.根据权利要求1所述基于球谐展开的GNSS标准单点定位方法，其特征在于，
所述步骤I.3.3中，基于最小二乘谱修正迭代公式计算X的解的过程如下：
基于最小二乘谱修正迭代解算X包括迭代①和迭代②两个迭代过程；其中，设定用于判断迭代①是否收敛的第一比较阈值，设定用于判断迭代②是否收敛的第二比较阈值；
迭代①第一次迭代时，设定K的初始值为1，将公式(11)得到的截断奇异值解
Figure FDA0002979540060000051
作为公式(14)第一次迭代时X的初值，并赋值给公式(14)等式右边的X；
每次迭代过程中，将上一次迭代时利用公式(14)得到的未知参数X的值，作为本次迭代过程中X的初值，并赋值给公式(14)等式右边的X；
基于公式(14)求解得到每次迭代过程中未知参数X的值；
将每次迭代过程中求解得到的X的值，与每次迭代时X的初值进行差值运算；
若经过差值运算，得到的差值大于第一比较阈值，则表示迭代不收敛，需要改正最小二乘谱修正迭代中的K值，即将K值递增1，继续上述迭代①过程；
若经过差值运算得到的差值小于或等于第一比较阈值，则跳出迭代①，进入迭代②；
迭代②将利用公式(14)求得的未知参数X的值与第二比较阈值进行比较；
若利用公式(14)求解得到的未知参数X的值大于第二比较阈值，则表示迭代②未收敛，此时将位置参数(dX,dY,dZ)分别加到测站的近似坐标(X0,Y0,Z0)上，更新测站的近似坐标，通过公式(5)对更新测站近似坐标后的标准单点定位误差方程进行线性化，通过公式(6)简化线性化后的标准单点定位误差方程，重新进入迭代①求解未知参数X的值；
若利用公式(14)求解得到的未知参数X的值小于或等于第二比较阈值，则表示迭代②收敛，跳出迭代；此时公式(14)求得的未知参数X的值，即未知参数X的解。
4.根据权利要求1所述基于球谐展开的GNSS标准单点定位方法，其特征在于，
所述步骤I.3.4中，基于最小二乘法求解未知参数X的过程如下：
设定用于判断迭代是否收敛的第三比较阈值；每次迭代过程中，利用公式(8)求解出未知参数X的值，将本次迭代得到的未知参数X的值与第三比较阈值进行比较；
若利用公式(8)求解得到的未知参数X的值大于第三比较阈值，则表示迭代未收敛；此时，将位置参数(dX,dY,dZ)分别加到测站的近似坐标(X0,Y0,Z0)上，更新测站的近似坐标，执行公式(5)-公式(8)的计算过程，并重新求解未知参数X的值；
若利用公式(8)求解得到的未知参数X的值小于或等于第三比较阈值，则表示迭代收敛，跳出迭代；此时，利用公式(8)得到的未知参数X的值，即待求解的未知参数X的解。
CN202110283670.9A 2021-03-17 2021-03-17 一种基于球谐展开的gnss单点定位方法 Active CN113093242B (zh)

Priority Applications (5)

Application Number Priority Date Filing Date Title
CN202110283670.9A CN113093242B (zh) 2021-03-17 2021-03-17 一种基于球谐展开的gnss单点定位方法
CN202210137675.5A CN114518586B (zh) 2021-03-17 2021-03-17 一种基于球谐展开的gnss精密单点定位方法
PCT/CN2021/123845 WO2022048694A1 (zh) 2021-03-17 2021-10-14 一种基于球谐展开的gnss单点定位方法
ZA2022/03722A ZA202203722B (en) 2021-03-17 2022-03-31 Gnss standard point positioning method based on spherical harmonics
US17/730,567 US20220299652A1 (en) 2021-03-17 2022-04-27 Gnss standard point positioning method based on spherical harmonics

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110283670.9A CN113093242B (zh) 2021-03-17 2021-03-17 一种基于球谐展开的gnss单点定位方法

Related Child Applications (1)

Application Number Title Priority Date Filing Date
CN202210137675.5A Division CN114518586B (zh) 2021-03-17 2021-03-17 一种基于球谐展开的gnss精密单点定位方法

Publications (2)

Publication Number Publication Date
CN113093242A CN113093242A (zh) 2021-07-09
CN113093242B true CN113093242B (zh) 2022-03-11

Family

ID=76668228

Family Applications (2)

Application Number Title Priority Date Filing Date
CN202210137675.5A Active CN114518586B (zh) 2021-03-17 2021-03-17 一种基于球谐展开的gnss精密单点定位方法
CN202110283670.9A Active CN113093242B (zh) 2021-03-17 2021-03-17 一种基于球谐展开的gnss单点定位方法

Family Applications Before (1)

Application Number Title Priority Date Filing Date
CN202210137675.5A Active CN114518586B (zh) 2021-03-17 2021-03-17 一种基于球谐展开的gnss精密单点定位方法

Country Status (4)

Country Link
US (1) US20220299652A1 (zh)
CN (2) CN114518586B (zh)
WO (1) WO2022048694A1 (zh)
ZA (1) ZA202203722B (zh)

Families Citing this family (18)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114518586B (zh) * 2021-03-17 2024-04-30 山东科技大学 一种基于球谐展开的gnss精密单点定位方法
CN114779285A (zh) * 2022-04-18 2022-07-22 浙江大学 基于微小型低功耗星载双模四频gnss接收机的精密定轨方法
CN114935768B (zh) * 2022-07-13 2022-11-04 武汉大学 一种基于单基站的虚拟参考站的构建方法
CN115168788B (zh) * 2022-09-07 2022-11-22 中国科学院空天信息创新研究院 卫星遥感大数据的确定方法、装置、设备及介质
CN115598673B (zh) * 2022-09-29 2023-10-24 同济大学 Igs gnss卫星钟差和轨道单天相邻产品边界处偏差计算方法
CN115616625B (zh) * 2022-10-08 2023-07-28 国家基础地理信息中心 一种gnss实时数据偏移方法及系统
CN116027459B (zh) * 2022-12-30 2024-05-14 桂林理工大学 基于数值天气预报数据的大气加权平均温度的计算方法
CN116029177B (zh) * 2023-03-14 2023-06-09 中国科学院空天信息创新研究院 一种融合球谐函数和精细网格的对流层延迟建模方法
BE1029960B1 (de) * 2023-04-06 2024-05-06 Univ Shandong Science & Tech GNSS-Standard-Einzelpunkt-Positionierungssystem und -verfahren
CN116108328B (zh) * 2023-04-13 2023-09-05 中国人民解放军63921部队 并置站不同天线参考点的相对位置获取方法和存储介质
CN116256788B (zh) * 2023-05-11 2023-07-11 中国人民解放军战略支援部队航天工程大学 一种基于阿波罗尼斯圆的空间几何迭代卫星定位方法
CN116609810B (zh) * 2023-05-19 2024-06-07 复旦大学 基于导航地基系统的电离层四维电子密度动态预测方法
CN116611329B (zh) * 2023-05-19 2024-05-03 复旦大学 基于深度算子网络的电离层电子总含量的四维估计方法
CN117421525B (zh) * 2023-12-18 2024-03-08 湖南科技大学 一种病态问题解算精度的相对均方误差分析方法
CN117454679B (zh) * 2023-12-26 2024-05-03 中铁第一勘察设计院集团有限公司 基于线状分布测站的动态线状电离层建模方法及系统
CN117492043B (zh) * 2023-12-29 2024-05-03 天津云遥宇航科技有限公司 一种基于shapiro时延和周期性相对论效应改正方法
CN117761740B (zh) * 2024-02-22 2024-05-14 开普勒卫星科技(武汉)有限公司 多系统基准站接收机精密脱敏算法
CN118013768A (zh) * 2024-04-10 2024-05-10 中国科学院地质与地球物理研究所 一种行星岩石圈磁场模型系数的确定方法及装置

Citations (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2014065664A1 (en) * 2012-10-25 2014-05-01 Fugro N.V. Ppp-rtk method and system for gnss signal based position determination
CN105182388A (zh) * 2015-10-10 2015-12-23 安徽理工大学 一种快速收敛的精密单点定位方法
EP2995975A1 (en) * 2014-09-15 2016-03-16 Fugro N.V. Precise gnss positioning system with improved ambiguity estimation
EP2995972A1 (en) * 2014-09-15 2016-03-16 Fugro N.V. Integer ambiguity-fixed precise point positioning method and system
CN106407560A (zh) * 2016-09-19 2017-02-15 武汉大学 表征大气各向异性的对流层映射函数模型的构建方法
CN107632313A (zh) * 2017-09-13 2018-01-26 航天恒星科技有限公司 基于相关性的卫星导航信号和sbas电文仿真方法
CN107765275A (zh) * 2017-09-04 2018-03-06 深圳市时空导航科技有限公司 广域差分定位方法、装置、终端及计算机可读存储介质
CN110275185A (zh) * 2019-07-11 2019-09-24 武汉大学 基于gnss和geo卫星的电离层投影函数建模方法
CN111551975A (zh) * 2020-06-24 2020-08-18 辽宁工程技术大学 Bds/gps参考站低高度角卫星整周模糊度确定方法
CN111551971A (zh) * 2020-05-14 2020-08-18 中国北方工业有限公司 一种支持异频gnss信号伪距差分定位的方法
CN111812681A (zh) * 2020-08-24 2020-10-23 中国人民解放军海军工程大学 大气区域建模方法、装置、电子设备及存储介质
CN112034500A (zh) * 2020-08-20 2020-12-04 上海华测导航技术股份有限公司 基于实时ppp模糊度固定技术的区域格网电离层建模方法

Family Cites Families (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102096084B (zh) * 2010-12-09 2012-12-26 东南大学 基于星间组合差分的精密单点定位方法
US9557419B2 (en) * 2012-12-18 2017-01-31 Trimble Inc. Methods for generating accuracy information on an ionosphere model for satellite navigation applications
CN103323888B (zh) * 2013-04-24 2015-06-17 东南大学 Gnss大气探测数据中对流层延迟误差的消除方法
EP3035080B1 (en) * 2014-12-16 2022-08-31 Trimble Inc. Navigation satellite system positioning involving the generation of correction information
EP3130943B1 (en) * 2015-08-14 2022-03-09 Trimble Inc. Navigation satellite system positioning involving the generation of tropospheric correction information
CN107356947B (zh) * 2017-05-31 2019-06-18 中国科学院测量与地球物理研究所 基于单频导航卫星数据确定卫星差分伪距偏差的方法
CN107153209B (zh) * 2017-07-06 2019-07-30 武汉大学 一种短弧段低轨导航卫星实时精密定轨方法
CN108845340A (zh) * 2018-06-01 2018-11-20 浙江亚特电器有限公司 基于gnss-rtk的定位方法
CN109709591B (zh) * 2018-12-07 2021-04-20 中国科学院光电研究院 一种面向智能终端的gnss高精度定位方法
CN109613582B (zh) * 2018-12-17 2021-11-23 中国科学院国家授时中心 一种车载实时单频米级伪距定位方法
CN114518586B (zh) * 2021-03-17 2024-04-30 山东科技大学 一种基于球谐展开的gnss精密单点定位方法

Patent Citations (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2014065664A1 (en) * 2012-10-25 2014-05-01 Fugro N.V. Ppp-rtk method and system for gnss signal based position determination
EP2995975A1 (en) * 2014-09-15 2016-03-16 Fugro N.V. Precise gnss positioning system with improved ambiguity estimation
EP2995972A1 (en) * 2014-09-15 2016-03-16 Fugro N.V. Integer ambiguity-fixed precise point positioning method and system
CN105182388A (zh) * 2015-10-10 2015-12-23 安徽理工大学 一种快速收敛的精密单点定位方法
CN106407560A (zh) * 2016-09-19 2017-02-15 武汉大学 表征大气各向异性的对流层映射函数模型的构建方法
CN107765275A (zh) * 2017-09-04 2018-03-06 深圳市时空导航科技有限公司 广域差分定位方法、装置、终端及计算机可读存储介质
CN107632313A (zh) * 2017-09-13 2018-01-26 航天恒星科技有限公司 基于相关性的卫星导航信号和sbas电文仿真方法
CN110275185A (zh) * 2019-07-11 2019-09-24 武汉大学 基于gnss和geo卫星的电离层投影函数建模方法
CN111551971A (zh) * 2020-05-14 2020-08-18 中国北方工业有限公司 一种支持异频gnss信号伪距差分定位的方法
CN111551975A (zh) * 2020-06-24 2020-08-18 辽宁工程技术大学 Bds/gps参考站低高度角卫星整周模糊度确定方法
CN112034500A (zh) * 2020-08-20 2020-12-04 上海华测导航技术股份有限公司 基于实时ppp模糊度固定技术的区域格网电离层建模方法
CN111812681A (zh) * 2020-08-24 2020-10-23 中国人民解放军海军工程大学 大气区域建模方法、装置、电子设备及存储介质

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
《Improvement of Indian-Regional Klobuchar Ionospheric Model Parameters for Single-Frequency GNSS Users》;D.Venkata Ratnam et.al;《IEEE GEOSCIENCE AND REMOTE SENSING LETTERS》;20180731;第15卷(第7期);第971-975页 *
《基于多通道奇异谱分析的GNSS坐标时间序列共模误差的提取》;周茂盛 等;《地球物理学报》;20181130;第61卷(第11期);第4383-4395页 *

Also Published As

Publication number Publication date
ZA202203722B (en) 2022-06-29
WO2022048694A1 (zh) 2022-03-10
CN113093242A (zh) 2021-07-09
CN114518586A (zh) 2022-05-20
CN114518586B (zh) 2024-04-30
US20220299652A1 (en) 2022-09-22

Similar Documents

Publication Publication Date Title
CN113093242B (zh) 一种基于球谐展开的gnss单点定位方法
Chen et al. Voxel-optimized regional water vapor tomography and comparison with radiosonde and numerical weather model
Luthcke et al. The 1-centimeter orbit: Jason-1 precision orbit determination using GPS, SLR, DORIS, and altimeter data special issue: Jason-1 calibration/validation
RU2479855C2 (ru) Зависящее от расстояния уменьшение ошибки при определении местоположения в режиме кинематики реального времени
CN108120994B (zh) 一种基于星载gnss的geo卫星实时定轨方法
Yao et al. Global ionospheric modeling based on multi-GNSS, satellite altimetry, and Formosat-3/COSMIC data
Paziewski Study on desirable ionospheric corrections accuracy for network-RTK positioning and its impact on time-to-fix and probability of successful single-epoch ambiguity resolution
CN112129300B (zh) 位置间动力学约束的低轨卫星星载gnss精密定轨方法及系统
CN107966722B (zh) 一种gnss钟差解算方法
Xu et al. Autonomous broadcast ephemeris improvement for GNSS using inter-satellite ranging measurements
CN115047499B (zh) 星载gnss-r土壤温度与湿度的反演方法、系统
Zhao et al. Accuracy and reliability of tropospheric wet refractivity tomography with GPS, BDS, and GLONASS observations
CN113902645A (zh) 一种基于逆rd定位模型的星载sar图像rpc校正参数获取方法
Li et al. Real‐time sensing of precipitable water vapor from beidou observations: hong kong and CMONOC networks
CN115755115A (zh) 基于gnss对流层层析技术的ppp改善方法
Zhou et al. Real-time orbit determination of Low Earth orbit satellite based on RINEX/DORIS 3.0 phase data and spaceborne GPS data
Brack et al. Operational multi-GNSS global ionosphere maps at GFZ derived from uncombined code and phase observations
Raquet et al. Use of a Covariance Analysis Technique for Predicting Performance of Regional‐Area Differential Code and Carrier‐Phase Networks
Martin GNSS precise point positioning: The enhancement with GLONASS
CN100371731C (zh) Gps和伪卫星组合定位方法
CN116594046B (zh) 基于低轨卫星信号多普勒误差补偿的运动目标定位方法
CN104991265A (zh) 一种北斗卫星导航系统用户统一性定位方法
CN115308781B (zh) 基于bdgim辅助的相位平滑伪距高精度时间传递方法
Li et al. A grid-based ionospheric weighted method for PPP-RTK with diverse network scales and ionospheric activity levels
Hong et al. Analysis of dual-frequency solution method for single-frequency precise point positioning based on SEID model for GPS and BDS

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant