CN116882134A - 一种多基准站网的gnss基线联合解算方法及计算机可读介质 - Google Patents
一种多基准站网的gnss基线联合解算方法及计算机可读介质 Download PDFInfo
- Publication number
- CN116882134A CN116882134A CN202310707640.5A CN202310707640A CN116882134A CN 116882134 A CN116882134 A CN 116882134A CN 202310707640 A CN202310707640 A CN 202310707640A CN 116882134 A CN116882134 A CN 116882134A
- Authority
- CN
- China
- Prior art keywords
- reference station
- double
- baseline
- difference
- ambiguity
- 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.)
- Pending
Links
- 238000004364 calculation method Methods 0.000 title claims description 15
- 238000000034 method Methods 0.000 claims abstract description 52
- 239000005436 troposphere Substances 0.000 claims abstract description 22
- 230000008569 process Effects 0.000 claims abstract description 15
- 238000001914 filtration Methods 0.000 claims abstract description 10
- 238000011084 recovery Methods 0.000 claims abstract description 5
- 230000008859 change Effects 0.000 claims abstract description 4
- 238000005295 random walk Methods 0.000 claims abstract description 4
- 239000011159 matrix material Substances 0.000 claims description 80
- 239000013598 vector Substances 0.000 claims description 46
- 238000013178 mathematical model Methods 0.000 claims description 22
- 238000012937 correction Methods 0.000 claims description 18
- 238000004422 calculation algorithm Methods 0.000 claims description 16
- 238000007781 pre-processing Methods 0.000 claims description 9
- 238000004590 computer program Methods 0.000 claims description 7
- 238000012545 processing Methods 0.000 claims description 7
- 238000001514 detection method Methods 0.000 claims description 6
- 239000000470 constituent Substances 0.000 claims description 3
- 230000008030 elimination Effects 0.000 claims description 3
- 238000003379 elimination reaction Methods 0.000 claims description 3
- 238000005259 measurement Methods 0.000 claims description 3
- 230000001360 synchronised effect Effects 0.000 abstract 1
- 238000012544 monitoring process Methods 0.000 description 7
- 238000005516 engineering process Methods 0.000 description 4
- 238000011160 research Methods 0.000 description 4
- 230000003068 static effect Effects 0.000 description 3
- 238000010586 diagram Methods 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000012805 post-processing Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
-
- 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/35—Constructional details or hardware or software details of the signal processing chain
- G01S19/37—Hardware or software details of the signal processing chain
-
- 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
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- Theoretical Computer Science (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Computational Mathematics (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Pure & Applied Mathematics (AREA)
- Data Mining & Analysis (AREA)
- General Engineering & Computer Science (AREA)
- Signal Processing (AREA)
- Computer Networks & Wireless Communication (AREA)
- Algebra (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- Operations Research (AREA)
- Computing Systems (AREA)
- Geometry (AREA)
- Evolutionary Computation (AREA)
- Computer Hardware Design (AREA)
- Position Fixing By Use Of Radio Waves (AREA)
Abstract
本发明公开了一种多基准站网的GNSS基线联合解算方法及计算机可读介质,本发明方法利用Delaunay三角形构建最优基线网,利用载波相位观测值对各基线逐个进行双差模糊度解算,恢复出各卫星的整周模糊度;利用所有基线中模糊度恢复后的载波相位观测值构建观测方程,方程中待估参数为各基准站坐标和天顶对流层湿延迟;同时构建准静态卡尔曼滤波模型,采用随机游走过程描述对流层延迟的状态变化;最后通过参数估计确定各基准站坐标的最优估值。提出的多基线联合解算方法,既充分考虑了同步基线之间的误差相关性,确保了解算方法在理论上的严密性,提高了平差解算结果的精度和可靠性,又大量减少了整体平差解算时模型参数的数量,提高了计算速度。
Description
技术领域
本发明涉及高精度卫星导航定位数据处理技术,尤其涉及一种多基准站网的GNSS基线联合解算方法及计算机可读介质。
背景技术
GNSS技术由于具有自动化、实时性、高精度和全天候等优点,成为地表变形监测的常用方法。对于大坝、桥梁等变形运动较为缓慢的监测场景来说,通常采用多静态基线网监测法。静态变形监测主要是通过对形变体进行周期性重复观测,掌握其变形规律,从而为地质灾害安全监测预警提供数据源。
在GNSS静态网监测中,基线解算技术作为高精度数据后处理的重要组成部分,是保障监测结果达到毫米级精度的关键。基线解算技术是多个测站共同观测多颗卫星利用测站之间的相关性,采用载波相位观测值的不同线性组合,解算测站之间的相对矢量。在基线解算相关方法的研究成果方面,周乐韬等对不同基线解算模式进行了分析,将多基线组成的星形网作为基本解算单元进行模糊度解算;唐卫明等对北斗导航系统单历元进行模糊度固定和基线解算研究,并进行精度对比分析;王振方等对GPS多基线模式下观测值的相关性进行了详细的分析;储立新在GPS/GLONASS组合的多基线解算中引入了等价方程的概念,用以组建多基线解算方程;祝会忠等研究了一种BDS中长距离基线高精度静态定位方法,在解算出宽巷整周模糊度的基础上,再对天顶对流层延迟误差进行参数估计,同时进行载波相位整周模糊度解算和定位计算。
上述研究工作在基线解算过程中主要采用单基线解算或多基线整体解算的模式,但是两种模式各自存在一定的局限性。单基线模式在独立基线解算的过程中,未能充分利用站点之间的相关性,理论上不够严密。而多基线整体解模式可以顾及参数间的关联性,数学模型更加严密,但是当站点数较多的情况下,解算参数比较多,数学模型和解算过程都比较复杂,在一定程度上限制了该方法的使用。
发明内容
为了解决上述技术问题,本发明提出了一种多基准站网的GNSS基线联合解算方法及计算机可读介质。
本发明方法的技术方案为一种多基准站网的GNSS基线联合解算方法,包括以下步骤:
步骤1:获取每个基准站的初始三维坐标,将每个基准站的初始三维坐标通过高斯投影算法转换得到每个基准站在高斯平面中的坐标,通过多个基准站在高斯平面中的坐标构建基准站平面坐标点集,将基准站平面坐标点集通过三角网生长算法得到多条待解算的基线;
步骤2:解算每条待解算基线的双差整周模糊度;
步骤3:依次构建每条待平差基线的误差方程、多基线三角网的误差方程、多基线三角网双差观测值权阵、多个基准站的平差数学模型,结合多个基准站的平差数学模型计算每个基准站的校正后三维坐标。
作为优选,步骤1所述基准站平面坐标点集,定义如下:
P={Pi|i=1,2,...,k}
其中,Pi表示第i个基准站在高斯平面中的坐标,k表示基准站的数量;
步骤1所述将基准站平面坐标点集通过三角网生长算法得到多条待解算的基线,具体步骤如下:
步骤1.1:在点集P中任意选择一点作为初始三角形的第一个顶点;
步骤1.2:搜索点集P中距离第一个顶点最近的点,作为初始三角形的第二个顶点,两点相连生成初始边;
步骤1.3:在点集P中搜索距离初始边中点最近且不和初始三角形的第一个顶点、初始三角形的第二个顶点在一条直线上的点作为初始三角形的第三个顶点,三个顶点连接成三角形,即为初始三角形;
步骤1.4:通过步骤1.3确定的两条边为基础,重复步骤1.3直到所有边构建完毕,此时所有边即可组成一个平面三角网;
则平面三角网中每条边对应两个基准站在三维空间中构成的向量即为每条待解算的基线。
作为优选,步骤2所述解算每条待解算基线的双差整周模糊度,具体如下:
步骤2.1:将每条待解算基线的GNSS载波相位观测进行数据预处理,得到预处理后两个基准站的GNSS载波相位观测数据;
选取每条待解算基线,假设基线两端的基准站分别为第i个基准站和第j个基准站,对这两个基准站的GNSS载波相位观测数据分别进行粗差探测和周跳探测,将含有粗差的GNSS载波相位观测数据进行粗差剔除处理,将含有周跳的GNSS载波相位观测数据进行周跳修复处理,得到组成该基线预处理后两个基准站的GNSS载波相位观测数据;
步骤2.2:利用误差传播定律确定双差观测值的权矩阵;
使用预处理后第i个基准站和第j个基准站的GNSS载波相位观测数据,构建每条待解算的基线经过线性化的载波相位双差观测方程;并在卫星高度角定权模型确定原始观测值权矩阵的基础上,再利用误差传播定律确定双差观测值的权矩阵P;
步骤2.3:解算每条待解算基线的双差整周模糊度;
将基准站i和j所组成基线经过线性化的载波相位双差观测方程转为误差方程,表示如下:
V=Aa+Bb-l
式中,V为双差观测值残差向量,A为对角线元素为λ的对角矩阵,a为双差整周模糊度向量;B为包含方向余弦及双差对流层投影函数在内的系数矩阵,b为由基准站j的坐标改正数dXj和天顶方向上两基准站间相对对流层湿延迟dTw构成的参数向量;l为载波相位双差观测值与其对应计算值之间的差值;
依据上述误差方程,忽略双差整周模糊度的整数特性,采用最小二乘平差求解双差模糊度向量估值基线及对流层参数向量/>及其协因数矩阵Q,其公式表示如下:
式中,为双差模糊度的协因数矩阵,/>为基线及对流层参数的协因数矩阵,与/>互为转置矩阵且均为模糊度参数与基线及对流层参数间的协因数矩阵;
利用双差模糊度向量估值及其协因数矩阵/>并结合最小二乘模糊度降相关平差算法,通过搜索得到最优的双差整周模糊度向量/>次优模糊度与最优模糊度残差二次型的比值R和模糊度解算成功率Ps;
若比值R和成功率Ps均大于设定阈值,则最优的双差整周模糊度向量即为所选取基线最终的双差模糊度固定值,此时该基线的双差整周模糊度解算完成;
重复执行步骤2.1-步骤2.3,直到所有待解算基线的双差整周模糊度解算完成。
作为优选,步骤3所述构建每条待平差基线的误差方程,具体如下:
将每条待解算基线视为每条待平差基线;
选取平面三角网中由第i个基准站和第j个基准站组成的待平差基线,将对应的模糊度固定值重新代入载波相位双差观测方程,并对构成基线的两个测站坐标进行线性化,经过整理可得到如下误差方程:
i,j∈{1,2,...,k}且i≠j
式中,为模糊度恢复后第i个基准站和第j个基准站之间的双差观测值残差向量,系数矩阵/>为第j个基准站中包含方向余弦及对流层投影函数在内的系数矩阵,/>为第i个基准站中包含方向余弦及对流层投影函数在内的系数矩阵,/>为第j个基准站的坐标改正数及天顶方向上对流层湿延迟dTw,j构成的参数向量,/>为第i个基准站的坐标改正数及天顶方向上对流层湿延迟dTw,i构成的参数向量,常数项/>为模糊度恢复后第i个基准站和第j个基准站之间双差观测值减去计算初值后的向量,k表示基准站个数。
步骤3所述构建多基线三角网的误差方程,具体如下:
叠加三角网中所有待平差基线的误差方程,构建多基线三角网观测值的误差方程,具体如下:
Vo=BoX-Lo
式中,Vo为所有待平差基线在模糊度恢复后的双差观测值的残差向量,Bo为所有待平差基线的误差方程叠加后待估参数X的系数矩阵,待估参数即由各基准站的坐标改正数及天顶方向上对流层湿延迟构成的参数向量/>组成,i=1,2,...,k,k表示基准站个数,常数项Lo为所有待平差基线在模糊度恢复后的双差观测值减去其计算初值后的向量。
步骤3所述构建多基线三角网双差观测值权阵,具体如下:
针对模糊度恢复后所有待平差基线的GNSS载波相位观测数据,采用卫星高度角定权模型确定原始观测值的方差矩阵,再利用误差传播定律确定双差观测值的方差-协方差矩阵,取逆阵即为观测值权矩阵Po。
步骤3所述构建坐标基准约束的基准站平差数学模型,具体如下:
若在三角网平差时需要引入坐标基准约束,则在步骤3.2的基础上增加约束方程:
Vc=TX
式中,Vc为坐标基准约束的残差,T为4c×4k的系数矩阵,c为约束基准的个数且c<k,k为基准站个数,矩阵T的定义如下:假设第m,m=1,2,...,k,个基准站为被坐标约束的基准站,其在所有被坐标约束的基准站中排序第n,n=1,2,...,c,则矩阵T中第(4n-3)行第(4m-3)列、第(4n-2)行第(4m-2)列及第(4n-1)行第(4m-1)列的元素为1,其余所有元素均为0;
在观测值权矩阵基础上,以追加分块对角矩阵的形式增加约束方程的权矩阵Pc,其维数为4c×4c;
矩阵Pc的定义如下:该矩阵中第(4n-3)、第(4n-2)及第(4n-1)个对角线元素,n=1,2,...,c,取极大值M,其他元素均取值0,则多个基准站的平差数学模型为:
步骤3所述结合多个基准站的平差数学模型计算每个基准站的校正后三维坐标,具体如下:
将多个基准站的平差数学模型作为卡尔曼滤波的观测方程,采用随机游走过程描述测站坐标和天顶对流层延迟参数的状态变化,构建待估参数的状态方程;通过卡尔曼滤波算法得到每个基准站的坐标改正数和每个基准站的对流层参数的滤波解,将每个基准站的初始三维坐标叠加每个基准站的坐标改正数,得到每个基准站的校正后三维坐标。
本发明还提供了一种计算机可读介质,所述计算机可读介质存储电子设备执行的计算机程序,当所述计算机程序在电子设备上运行时,执行所述多基准站网的GNSS基线联合解算方法的步骤。
本发明产生的有益效果是:
本发明方法在传统单基线解算获取测站坐标和模糊度参数的基础上,在网平差时重新构建三角网各基线的观测方程对测站坐标参数进行统一解算。通过顾及观测值及各参数间的相关性,增加了平差模型的强度,提高了平差解算结果的精度和可靠性;
本发明方法在单基线模式恢复出各卫星对整周模糊度的基础上,利用多基线整体解算模式进行坐标和天顶对流层参数的联合估计。通过第一步的模糊度解算,大量减少了基线网平差时待估参数的数量,提高了三角网平差解算的速度。
附图说明
图1:本发明实施例的点位分布示意图;
图2:本发明实施例的方法流程图。
具体实施方式
下面将结合本发明实施例中的附图,对本发明实施例中的技术方案进行清楚、完整地描述,显然,所描述的实施例仅仅是本发明一部分实施例,而不是全部的实施例。基于本发明中的实施例,本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例,都属于本发明保护的范围。
具体实施时,本发明技术方案提出的方法可由本领域技术人员采用计算机软件技术实现自动运行流程,实现方法的系统装置例如存储本发明技术方案相应计算机程序的计算机可读存储介质以及包括运行相应计算机程序的计算机设备,也应当在本发明的保护范围内。
为了使本发明的目的、技术方案及优点更加清楚明白,以下结合4个基准站的实施例,对本发明进行进一步详细说明。四个基准站的点位分布示意图如图1所示,其中1和3为坐标固定且已知的基准站,在后续平差处理中将当作约束基准。应当理解,此处所描述的具体实施例仅用以解释本发明,并不用于限定本发明。
下面结合图1-2介绍本发明实施例的技术方案为一种多基准站网的GNSS基线联合解算方法,具体如下:
如图2所示,为本发明实施例的方法流程图。
步骤1:获取每个基准站的初始三维坐标,将每个基准站的初始三维坐标通过高斯投影算法转换得到每个基准站在高斯平面中的坐标,通过多个基准站在高斯平面中的坐标构建基准站平面坐标点集,将基准站平面坐标点集通过三角网生长算法得到多条待解算的基线;
步骤2:解算每条待解算基线的双差整周模糊度;
步骤3:依次构建每条待平差基线的误差方程、多基线三角网的误差方程、多基线三角网双差观测值权阵、多个基准站的平差数学模型,结合多个基准站的平差数学模型计算每个基准站的校正后三维坐标;
步骤1所述基准站平面坐标点集,定义如下:
P={Pi|i=1,2,...,k}
其中,Pi表示第i个基准站在高斯平面中的坐标,k=4表示基准站的数量;
步骤1所述将基准站平面坐标点集通过三角网生长算法得到多条待解算的基线,具体步骤如下:
步骤1.1:在点集P中任意选择点1作为初始三角形的第一个顶点;
步骤1.2:搜索点集P中距离点1最近的点2,作为初始三角形的第二个顶点,两点相连生成初始边L12;
步骤1.3:在点集P中搜索距离初始边L12中点最近且不和点1、点2在一条直线上的点3作为初始三角形的第三个顶点,三个顶点连接成三角形,即为初始三角形;
步骤1.4:通过步骤1.3确定的两条边L23和L13为基础,重复步骤1.3直到所有边构建完毕,此时所有边即可组成一个平面三角网;
则平面三角网中每条边对应两个基准站在三维空间中构成的向量即为每条待解算的基线;
步骤2所述解算每条待解算基线的双差整周模糊度,具体如下:
步骤2.1:将每条待解算基线的GNSS载波相位观测进行数据预处理,得到预处理后两个基准站的GNSS载波相位观测数据;
选取每条待解算基线,假设基线两端的基准站分别为基准站1和基准站2,对这两个基准站的GNSS载波相位观测数据分别进行粗差探测和周跳探测,将含有粗差的GNSS载波相位观测数据进行粗差剔除处理,将含有周跳的GNSS载波相位观测数据进行周跳修复处理,得到组成该基线预处理后两个基准站的GNSS载波相位观测数据;
步骤2.2:利用误差传播定律确定双差观测值的权矩阵;
使用预处理后基准站1和基准站2的GNSS载波相位观测数据,构建每条待解算的基线经过线性化的载波相位双差观测方程;并在卫星高度角定权模型确定原始观测值权矩阵的基础上,再利用误差传播定律确定双差观测值的权矩阵P;
步骤2.3:解算每条待解算基线的双差整周模糊度;
将基准站1和基准站2所组成基线经过线性化的载波相位双差观测方程转为误差方程,表示如下:
V=Aa+Bb-l
式中,V为双差观测值残差向量,A为对角线元素为λ的对角矩阵,a为双差整周模糊度向量;B为包含方向余弦及双差对流层投影函数在内的系数矩阵,b为由基准站2的坐标改正数dX2和天顶方向上两基准站间相对对流层湿延迟dTw构成的参数向量;l为载波相位双差观测值与其对应计算值之间的差值;
依据上述误差方程,忽略双差整周模糊度的整数特性,采用最小二乘平差求解双差模糊度向量估值基线及对流层参数向量/>及其协因数矩阵Q,其公式表示如下:
式中,为双差模糊度的协因数矩阵,/>为基线及对流层参数的协因数矩阵,与/>互为转置矩阵且均为模糊度参数与基线及对流层参数间的协因数矩阵;
利用双差模糊度向量估值及其协因数矩阵/>并结合最小二乘模糊度降相关平差算法,通过搜索得到最优的双差整周模糊度向量/>次优模糊度与最优模糊度残差二次型的比值R和模糊度解算成功率Ps。若比值R和成功率Ps均大于设定阈值,则最优的双差整周模糊度向量即为所选取基线最终的双差模糊度固定值,此时该基线的双差整周模糊度解算完成;
重复执行步骤2.1-步骤2.3,直到所有待解算基线的双差整周模糊度解算完成;
步骤3所述构建每条待平差基线的误差方程,具体如下:
将每条待解算基线视为每条待平差基线;
选取平面三角网中由第i个基准站和第j个基准站组成的待平差基线,将对应的模糊度固定值重新代入载波相位双差观测方程,并对构成基线的两个测站坐标进行线性化,经过整理可得到如下误差方程:
i,j∈{1,2,...,k}且i≠j
式中,为模糊度恢复后第i个基准站和第j个基准站之间的双差观测值残差向量,系数矩阵/>为第j个基准站中包含方向余弦及对流层投影函数在内的系数矩阵,/>为第i个基准站中包含方向余弦及对流层投影函数在内的系数矩阵,/>为第j个基准站的坐标改正数及天顶方向上对流层湿延迟dTw,j构成的参数向量,/>为第i个基准站的坐标改正数及天顶方向上对流层湿延迟dTw,i构成的参数向量,常数项/>为模糊度恢复后第i个基准站和第j个基准站之间双差观测值减去计算初值后的向量,k=4表示基准站个数,对于本例(i,j)可取(1,2),(1,3),(1,4),(2,3),(3,4);
步骤3所述构建多基线三角网的误差方程,具体如下:
叠加三角网中所有待平差基线的误差方程,构建多基线三角网观测值的误差方程,具体如下:
Vo=BoX-Lo
式中,Vo为所有待平差基线在模糊度恢复后的双差观测值的残差向量,Bo为所有待平差基线的误差方程叠加后待估参数X的系数矩阵,待估参数即由各基准站的坐标改正数及天顶方向上对流层湿延迟构成的参数向量/>组成,i=1,2,...,k,k=4表示基准站个数,常数项Lo为所有待平差基线在模糊度恢复后的双差观测值减去其计算初值后的向量;对于本例,
步骤3所述构建多基线三角网双差观测值权阵,具体如下:
针对模糊度恢复后所有待平差基线的GNSS载波相位观测数据,采用卫星高度角定权模型确定原始观测值的方差矩阵,再利用误差传播定律确定双差观测值的方差-协方差矩阵,取逆阵即为观测值权矩阵Po;
步骤3所述构建多个基准站的平差数学模型,具体如下:
若在三角网平差时需要引入坐标基准约束,则在步骤3.2的基础上增加约束方程:
Vc=TX
式中,Vc为坐标基准约束的残差,T为4c×4k的系数矩阵,c=2为约束基准的个数,k=4为基准站个数,矩阵T的定义如下:假设第m,m=1,2,...,k,个基准站为被坐标约束的基准站,其在所有被坐标约束的基准站中排序第n,n=1,2,...,c,则矩阵T中第(4n-3)行第(4m-3)列、第(4n-2)行第(4m-2)列及第(4n-1)行第(4m-1)列的元素为1,其余所有元素均为0。对于本例,
式中,04×4表示一个4×4维的零值矩阵,表示一个4×4维的对角矩阵且前三个对角线元素为1,即
在观测值权矩阵基础上,以追加分块对角矩阵的形式增加约束方程的权矩阵Pc,其维数为4c×4c,c=2为约束基准的个数;
矩阵Pc的定义如下:该矩阵中第(4n-3)、第(4n-2)及第(4n-1)个对角线元素,n=1,2,...,c,取极大值M=100000,其他元素均取值0。对于本例:
则多个基准站的平差数学模型为:
步骤3所述结合多个基准站的平差数学模型计算每个基准站的校正后三维坐标,具体如下:
将多个基准站的平差数学模型作为卡尔曼滤波的观测方程,采用随机游走过程描述测站坐标和天顶对流层延迟参数的状态变化,构建待估参数的状态方程;通过卡尔曼滤波算法得到每个基准站的坐标改正数和每个基准站的对流层参数的滤波解,将每个基准站的初始三维坐标叠加每个基准站的坐标改正数,得到每个基准站的校正后三维坐标。
本发明的具体实施例还提供了一种计算机可读介质。
所述计算机可读介质为服务器工作站;
所述服务器工作站存储电子设备执行的计算机程序,当所述计算机程序在电子设备上运行时,使得所述电子设备执行本发明实施例的多基准站网的GNSS基线联合解算方法的步骤。
应当理解的是,本说明书未详细阐述的部分均属于现有技术。
应当理解的是,上述针对较佳实施例的描述较为详细,并不能因此而认为是对本发明专利保护范围的限制,本领域的普通技术人员在本发明的启示下,在不脱离本发明权利要求所保护的范围情况下,还可以做出替换或变形,均落入本发明的保护范围之内,本发明的请求保护范围应以所附权利要求为准。
Claims (10)
1.一种多基准站网的GNSS基线联合解算方法,其特征在于:
获取每个基准站的初始三维坐标,构建基准站平面坐标点集,通过三角网生长算法得到多条待解算的基线;
解算每条待解算基线的双差整周模糊度;
构建多个基准站的平差数学模型,计算每个基准站的校正后三维坐标。
2.根据权利要求1所述的多基准站网的GNSS基线联合解算方法,其特征在于,包括以下步骤:
步骤1:获取每个基准站的初始三维坐标,将每个基准站的初始三维坐标通过高斯投影算法转换得到每个基准站在高斯平面中的坐标,通过多个基准站在高斯平面中的坐标构建基准站平面坐标点集,将基准站平面坐标点集通过三角网生长算法得到多条待解算的基线;
步骤2:解算每条待解算基线的双差整周模糊度;
步骤3:依次构建每条待平差基线的误差方程、多基线三角网的误差方程、多基线三角网双差观测值权阵、多个基准站的平差数学模型,结合多个基准站的平差数学模型计算每个基准站的校正后三维坐标。
3.根据权利要求2所述的多基准站网的GNSS基线联合解算方法,其特征在于:
步骤1所述基准站平面坐标点集,定义如下:
P={Pi|i=1,2,…,k}
其中,Pi表示第i个基准站在高斯平面中的坐标,k表示基准站的数量;
步骤1所述将基准站平面坐标点集通过三角网生长算法得到多条待解算的基线,具体步骤如下:
步骤1.1:在点集P中任意选择一点作为初始三角形的第一个顶点;
步骤1.2:搜索点集P中距离第一个顶点最近的点,作为初始三角形的第二个顶点,两点相连生成初始边;
步骤1.3:在点集P中搜索距离初始边中点最近且不和初始三角形的第一个顶点、初始三角形的第二个顶点在一条直线上的点作为初始三角形的第三个顶点,三个顶点连接成三角形,即为初始三角形;
步骤1.4:通过步骤1.3确定的两条边为基础,重复步骤1.3直到所有边构建完毕,此时所有边即可组成一个平面三角网;
则平面三角网中每条边对应两个基准站在三维空间中构成的向量即为每条待解算的基线。
4.根据权利要求3所述的多基准站网的GNSS基线联合解算方法,其特征在于:
所述解算每条待解算基线的双差整周模糊度,具体如下:
步骤2.1:将每条待解算基线的GNSS载波相位观测进行数据预处理,得到预处理后两个基准站的GNSS载波相位观测数据;
选取每条待解算基线,假设基线两端的基准站分别为第i个基准站和第j个基准站,对这两个基准站的GNSS载波相位观测数据分别进行粗差探测和周跳探测,将含有粗差的GNSS载波相位观测数据进行粗差剔除处理,将含有周跳的GNSS载波相位观测数据进行周跳修复处理,得到组成该基线预处理后两个基准站的GNSS载波相位观测数据;
步骤2.2:利用误差传播定律确定双差观测值的权矩阵;
使用预处理后第i个基准站和第j个基准站的GNSS载波相位观测数据,构建每条待解算的基线经过线性化的载波相位双差观测方程;并在卫星高度角定权模型确定原始观测值权矩阵的基础上,再利用误差传播定律确定双差观测值的权矩阵P;
步骤2.3:解算每条待解算基线的双差整周模糊度;
将基准站i和j所组成基线经过线性化的载波相位双差观测方程转为误差方程,表示如下:
V=a+Bb-l
式中,V为双差观测值残差向量,A为对角线元素为λ的对角矩阵,a为双差整周模糊度向量;B为包含方向余弦及双差对流层投影函数在内的系数矩阵,b为由基准站j的坐标改正数dXj和天顶方向上两基准站间相对对流层湿延迟dTw构成的参数向量;l为载波相位双差观测值与其对应计算值之间的差值;
依据上述误差方程,忽略双差整周模糊度的整数特性,采用最小二乘平差求解双差模糊度向量估值基线及对流层参数向量/>及其协因数矩阵Q,其公式表示如下:
式中,为双差模糊度的协因数矩阵,/>为基线及对流层参数的协因数矩阵,/>与互为转置矩阵且均为模糊度参数与基线及对流层参数间的协因数矩阵;
利用双差模糊度向量估值及其协因数矩阵/>并结合最小二乘模糊度降相关平差算法,通过搜索得到最优的双差整周模糊度向量/>次优模糊度与最优模糊度残差二次型的比值R和模糊度解算成功率Ps;
若比值R和成功率Ps均大于设定阈值,则最优的双差整周模糊度向量即为所选取基线最终的双差模糊度固定值,此时该基线的双差整周模糊度解算完成;
重复执行步骤2.1-步骤2.3,直到所有待解算基线的双差整周模糊度解算完成。
5.根据权利要求4所述的多基准站网的GNSS基线联合解算方法,其特征在于:
步骤3所述构建每条待平差基线的误差方程,具体如下:
将每条待解算基线视为每条待平差基线;
选取平面三角网中由第i个基准站和第j个基准站组成的待平差基线,将对应的模糊度固定值重新代入载波相位双差观测方程,并对构成基线的两个测站坐标进行线性化,经过整理可得到如下误差方程:
i,j∈{1,2,…,k}且i≠j
式中,为模糊度恢复后第i个基准站和第j个基准站之间的双差观测值残差向量,系数矩阵/>为第j个基准站中包含方向余弦及对流层投影函数在内的系数矩阵,/>为第i个基准站中包含方向余弦及对流层投影函数在内的系数矩阵,/>为第j个基准站的坐标改正数及天顶方向上对流层湿延迟dTw,j构成的参数向量,/>为第i个基准站的坐标改正数及天顶方向上对流层湿延迟dTw,i构成的参数向量,常数项/>为模糊度恢复后第i个基准站和第j个基准站之间双差观测值减去计算初值后的向量,k表示基准站个数。
6.根据权利要求5所述的多基准站网的GNSS基线联合解算方法,其特征在于:
步骤3所述构建多基线三角网的误差方程,具体如下:
叠加三角网中所有待平差基线的误差方程,构建多基线三角网观测值的误差方程,具体如下:
Vo=oX-o
式中,Vo为所有待平差基线在模糊度恢复后的双差观测值的残差向量,Bo为所有待平差基线的误差方程叠加后待估参数X的系数矩阵,待估参数即由各基准站的坐标改正数及天顶方向上对流层湿延迟构成的参数向量/>组成,i=1,2,…,k,k表示基准站个数,常数项Lo为所有待平差基线在模糊度恢复后的双差观测值减去其计算初值后的向量。
7.根据权利要求6所述的多基准站网的GNSS基线联合解算方法,其特征在于:
步骤3所述构建多基线三角网双差观测值权阵,具体如下:
针对模糊度恢复后所有待平差基线的GNSS载波相位观测数据,采用卫星高度角定权模型确定原始观测值的方差矩阵,再利用误差传播定律确定双差观测值的方差-协方差矩阵,取逆阵即为观测值权矩阵Po。
8.根据权利要求7所述的多基准站网的GNSS基线联合解算方法,其特征在于:
步骤3所述构建坐标基准约束的基准站平差数学模型,具体如下:
若在三角网平差时需要引入坐标基准约束,则在步骤3.2的基础上增加约束方程:
Vc=X
式中,Vc为坐标基准约束的残差,T为4×4k的系数矩阵,c为约束基准的个数且c<k,k为基准站个数,矩阵T的定义如下:假设第m,m=1,2,…,k,个基准站为被坐标约束的基准站,其在所有被坐标约束的基准站中排序第n,n=1,2,…,c,则矩阵T中第(4n-3)行第(4m-3)列、第(4n-2)行第(4m-2)列及第(4n-1)行第(4m-1)列的元素为1,其余所有元素均为0;
在观测值权矩阵基础上,以追加分块对角矩阵的形式增加约束方程的权矩阵Pc,其维数为4×4c;
矩阵Pc的定义如下:该矩阵中第(4n-3)、第(4n-2)及第(4n-1)个对角线元素,n=1,2,…,c,取极大值M,其他元素均取值0,则多个基准站的平差数学模型为:
9.根据权利要求8所述的多基准站网的GNSS基线联合解算方法,其特征在于:
步骤3所述结合多个基准站的平差数学模型计算每个基准站的校正后三维坐标,具体如下:
将多个基准站的平差数学模型作为卡尔曼滤波的观测方程,采用随机游走过程描述测站坐标和天顶对流层延迟参数的状态变化,构建待估参数的状态方程;通过卡尔曼滤波算法得到每个基准站的坐标改正数和每个基准站的对流层参数的滤波解,将每个基准站的初始三维坐标叠加每个基准站的坐标改正数,得到每个基准站的校正后三维坐标。
10.一种计算机可读介质,其特征在于,其存储电子设备执行的计算机程序,当所述计算机程序在电子设备上运行时,使得所述电子设备执行如权利要求1-9任一项所述方法的步骤。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202310707640.5A CN116882134A (zh) | 2023-06-14 | 2023-06-14 | 一种多基准站网的gnss基线联合解算方法及计算机可读介质 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202310707640.5A CN116882134A (zh) | 2023-06-14 | 2023-06-14 | 一种多基准站网的gnss基线联合解算方法及计算机可读介质 |
Publications (1)
Publication Number | Publication Date |
---|---|
CN116882134A true CN116882134A (zh) | 2023-10-13 |
Family
ID=88261221
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202310707640.5A Pending CN116882134A (zh) | 2023-06-14 | 2023-06-14 | 一种多基准站网的gnss基线联合解算方法及计算机可读介质 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN116882134A (zh) |
Cited By (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN117232380A (zh) * | 2023-11-13 | 2023-12-15 | 武汉大学 | 一种带状工程变形监测方法、装置、设备及可读存储介质 |
-
2023
- 2023-06-14 CN CN202310707640.5A patent/CN116882134A/zh active Pending
Cited By (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN117232380A (zh) * | 2023-11-13 | 2023-12-15 | 武汉大学 | 一种带状工程变形监测方法、装置、设备及可读存储介质 |
CN117232380B (zh) * | 2023-11-13 | 2024-01-30 | 武汉大学 | 一种带状工程变形监测方法、装置、设备及可读存储介质 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN114518586B (zh) | 一种基于球谐展开的gnss精密单点定位方法 | |
CN102288978B (zh) | 一种cors基站周跳探测与修复方法 | |
CN109917356B (zh) | 一种机载激光扫描系统误差标定方法 | |
Kampes et al. | Ambiguity resolution for permanent scatterer interferometry | |
Janssen et al. | Tropospheric corrections to SAR interferometry from GPS observations | |
CN104714244A (zh) | 一种基于抗差自适应Kalman滤波的多系统动态PPP解算方法 | |
CN111650579B (zh) | 一种岩移参数自适应获取的InSAR矿区三维形变估计方法、装置及介质 | |
CN110646822B (zh) | 一种基于惯导辅助的整周模糊度Kalman滤波算法 | |
CN116882134A (zh) | 一种多基准站网的gnss基线联合解算方法及计算机可读介质 | |
Sandwell et al. | Topographic phase recovery from stacked ERS interferometry and a low‐resolution digital elevation model | |
Wang et al. | Performance evaluation of a real-time high-precision landslide displacement detection algorithm based on GNSS virtual reference station technology | |
CN115877421A (zh) | 一种输电通道地质敏感区的形变检测方法及装置 | |
CN103760586B (zh) | 一种在gps姿态测量中快速探测与修复周跳的方法 | |
CN103760582B (zh) | 一种遮挡环境下卫星双差观测结构的优化方法 | |
CN110632636B (zh) | 一种基于Elman神经网络的载体姿态估计方法 | |
Ji et al. | Applying InSAR and GNSS data to obtain 3-D surface deformations based on iterated almost unbiased estimation and Laplacian smoothness constraint | |
Estahbanati et al. | A phase unwrapping approach based on extended Kalman filter for subsidence monitoring using persistent scatterer time series interferometry | |
Kuang et al. | Robust constrained Kalman filter algorithm considering time registration for GNSS/acoustic joint positioning | |
Xu et al. | A robust method for 3-D surface displacement fields combining GNSS and single-orbit InSAR measurements with directional constraint from elasticity model | |
Jokinen et al. | Improving fixed-ambiguity Precise Point Positioning (PPP) convergence time and accuracy by using GLONASS | |
Bozsó et al. | Integration of Sentinel-1 interferometry and GNSS networks for derivation of 3-D surface changes | |
King | The GPS contribution to the error budget of surface elevations derived from airborne LIDAR | |
Alves et al. | Introduction of a geometry-based Network RTK quality indicator | |
Awange | Diagnosis of Outlier of type multipath in GPS Pseudo-range observations | |
CN115390117B (zh) | 一种精密单点定位变形监测方法及装置 |
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 |