CN116204756B - 一种多分析中心精密站坐标产品综合方法及系统 - Google Patents

一种多分析中心精密站坐标产品综合方法及系统 Download PDF

Info

Publication number
CN116204756B
CN116204756B CN202310480488.1A CN202310480488A CN116204756B CN 116204756 B CN116204756 B CN 116204756B CN 202310480488 A CN202310480488 A CN 202310480488A CN 116204756 B CN116204756 B CN 116204756B
Authority
CN
China
Prior art keywords
station
coordinate
analysis
station coordinate
analysis center
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
CN202310480488.1A
Other languages
English (en)
Other versions
CN116204756A (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.)
Wuhan University WHU
Original Assignee
Wuhan University WHU
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 Wuhan University WHU filed Critical Wuhan University WHU
Priority to CN202310480488.1A priority Critical patent/CN116204756B/zh
Publication of CN116204756A publication Critical patent/CN116204756A/zh
Application granted granted Critical
Publication of CN116204756B publication Critical patent/CN116204756B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • 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
    • 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/396—Determining accuracy or reliability of position or pseudorange measurements
    • 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/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • 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/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Abstract

本发明公开了一种多分析中心精密站坐标产品综合方法及系统，所述方法包括：对多个分析中心提供的GNSS站坐标产品文件进行预处理，对原始法方程系统进行初步先验约束消除，利用数学变换消除残余约束；对目标框架进行非线性改正，将分析中心站坐标产品对齐到目标框架；建立站坐标综合模型，利用最小二乘原理得到综合解估计值和相似变换参数以及站坐标残差；识别分析中心站坐标产品中的异常站点，并进行各个分析中心产品的权重分配；迭代运算直至各个分析中心的权重趋于稳定或达到预设的最大迭代次数，求出综合站坐标。本发明综合利用了先验残余约束消除策略、框架对齐策略和基于站坐标相关性的粗差探测，具有高可靠性、高稳定性和高精度的优势。

Description

一种多分析中心精密站坐标产品综合方法及系统
技术领域
本发明属于全球卫星导航系统（GNSS）高精度定位技术领域，具体涉及基于iGMAS（全球连续监测评估系统）/IGS（国际GNSS服务组织）的一种多分析中心精密站坐标产品综合方法及系统。
背景技术
相对于单分析中心解算的站坐标产品，基于多分析中心站坐标产品的综合结果能极大改善站点坐标解的精度和可靠性，改善用户的产品使用体验。从成立以来，IGS陆续开展了多分析中心的站坐标产品的综合处理工作，中国正在建设和发展的全球连续监测评估系统（iGMAS）由30个跟踪站、3个数据中心、10余个分析中心、1个产品综合与服务中心和监测评估中心等组成，主要对GPS、GLONASS、BDS、Galileo四系统状态进行监测评估，并提供四系统高精度综合产品服务。高精度的综合产品不仅可以用来作为卫星状态监测评估的参考产品，评估不同卫星系统的广播轨道精度；还可以作为地球参考框架的二级实现，为地面站精密坐标获取提供基准，为研究地球科学提供重要的数据产品基础。
目前已有多个国内外机构提供全球GNSS跟踪站的坐标产品，也出现了一些多分析中心站坐标产品的综合方法，比如陈国等在《多分析中心站坐标产品的综合方法研究》中提出了站坐标和地球自转参数同时综合的方法，但是多分析中心站坐标产品的综合仍然没有较好的解决如下两个技术难点或问题：
其一，如何确保综合的站坐标产品最大程度上不受异常站点的影响。由于分析中心产品的质量受到观测数据质量及数据质量控制策略等影响，不同分析中心的站坐标产品在质量上存在差异，甚至可能会出现异常的站坐标解算结果，如何在站坐标产品综合中识别异常分析中心或者异常站点坐标，是提高站坐标综合结果可靠性的关键之一。
其二，如何将综合的站坐标产品对齐到目标框架如国际地球参考框架。目前，由于不同个分析中心在具体站坐标产品生成策略上（如模糊度固定策略、非差和双差观测量、先验约束等）存在差异，造成不同分析中心站坐标产品间的框架不一致，如何将相应的综合站坐标结果对齐到国际地球参考框架，是利用综合站坐标产品进行长时间序列分析的前提。上述两个技术难点直接影响综合站坐标产品的精度、稳定性和基准连续性，进而影响用户对综合产品的使用体验，是iGMAS/IGS的高精度产品服务中必须要解决的问题。
发明内容
有鉴于此，本发明提出了一种多分析中心精密站坐标产品综合方法及系统，用于解决现有的iGMAS/IGS产品精度难以保证的问题。
本发明第一方面，公开一种多分析中心精密站坐标产品综合方法，所述方法包括：
S1、对多个分析中心提供的GNSS站坐标产品文件进行预处理，利用分析中心提供的先验信息对预处理得到的原始法方程系统进行初步先验约束消除，利用相似变换矩阵消除残余约束；
S2、建立相似变换方程，采用非潮汐海洋、大气和水文负载产品对目标框架进行非线性改正，将消除先验约束之后的分析中心站坐标产品对齐到目标框架，并计算目标框架下的分析中心站坐标解；
S3、对于所有的分析中心站坐标产品，根据分析中心站坐标解建立站坐标综合模型，得到总的综合法方程系统，利用最小二乘原理得到综合解估计值和相似变换参数以及站坐标残差，实现多分析中心站坐标产品综合；
S4、利用站点间的相关性信息识别分析中心站坐标产品中的异常站点，并利用方差分量估计方法进行各个分析中心产品的权重分配；
S5、重复以上步骤S3~S4，直至各个分析中心的权重趋于稳定或达到预设的最大迭代次数，得到权重分配结果；
S6、基于权重分配结果求解总的综合法方程系统，输出综合站坐标。
在以上技术方案的基础上，优选的，所述预处理具体包括：初步判定各精密站坐标产品的格式，并统计精密站坐标产品中测站数、核心站数以及存储法方程系统信息。
在以上技术方案的基础上，优选的，所述利用相似变换矩阵消除残余约束具体包括：
设添加先验约束的站坐标估计系统如下：
Figure SMS_1
其中，
Figure SMS_2
为站坐标先验值，/>
Figure SMS_3
为站坐标估计值，/>
Figure SMS_4
为先验约束矩阵，/>
Figure SMS_5
为后验方差-协方差阵，/>
Figure SMS_6
为未添加任何约束的信息矩阵，u为原始法方程系统对应的右向量；
初步先验约束消除采用的公式如下：
Figure SMS_7
其中，
Figure SMS_8
为没有完全消除先验约束的信息矩阵；
采用下式进行残余约束的消除：
Figure SMS_9
其中，
Figure SMS_10
为单位矩阵，/>
Figure SMS_11
为无约束法方程系统的右矩阵，/>
Figure SMS_12
为无约束对应的坐标估计值；N为原始法方程系统中的信息矩阵，/>
Figure SMS_13
表示有n个数值1组成的列向量，n为跟踪站数，/>
Figure SMS_14
为克罗内克积，M为相似变换矩阵，M的表达式为：
Figure SMS_15
其中x i 、y i 、z i 分别为测站i的X、Y和Z坐标分量。
在以上技术方案的基础上，优选的，所述步骤S2具体包括如下分步骤：
计算站点在当前历元t的非潮汐海洋、大气和水文负载形变，非潮汐海洋、大气和水文负载形变记为
Figure SMS_16
；
在目标框架中引入非潮汐海洋、大气和水文负载，计算目标框架在当前历元的坐标，计算公式为：
Figure SMS_17
t和t 0分别表示当前历元和参考历元；
Figure SMS_18
表示目标框架解在当前历元的位置；
Figure SMS_19
表示目标框架解在参考历元的位置；/>
Figure SMS_20
表示目标框架解在参考历元的速度；
计算分析中心站坐标产品中核心站点目标框架解的坐标值
Figure SMS_21
；
对分析中心无约束的法方程系统添加统一约束：
Figure SMS_22
其中，
Figure SMS_23
，是由核心站坐标构成的约束矩阵，
Figure SMS_24
为核心站坐标构成的相似变换矩阵，/>
Figure SMS_25
为相应的转置矩阵，/>
Figure SMS_26
为目标框架下的分析中心站坐标解，s表示分析中心编号；
计算目标框架下的分析中心站坐标解
Figure SMS_27
。
在以上技术方案的基础上，优选的，所述站坐标综合模型的公式为：
Figure SMS_28
其中，s表示分心中心编号，
Figure SMS_29
表示分析中心站坐标残差，/>
Figure SMS_30
表示综合站坐标解；/>
Figure SMS_31
表示相似变换中的平移参数向量；/>
Figure SMS_32
表示相似变换中的旋转参数矩阵；/>
Figure SMS_33
相似变换中的尺度参数标量。
在以上技术方案的基础上，优选的，所述步骤S4中，所述利用站点间的相关性信息识别分析中心站坐标产品中的异常站点具体包括：
将每个分析中心的站坐标产品划分为一类观测数据，采用下式计算每个分析中心的后验方差：
Figure SMS_34
其中，
Figure SMS_35
表示后验方差；nsit表示分析中心s的站坐标产品中测站总数；设
Figure SMS_36
表示分析中心s中的跟踪站isit的坐标后验残差，/>
Figure SMS_37
表示分析中心s中的跟踪站k的坐标后验残差，/>
Figure SMS_38
表示站点isit和k之间的相关信息；/>
Figure SMS_39
表示综合法方程系统的信息矩阵逆；/>
Figure SMS_40
表示分析中心站坐标观测量的权阵；/>
Figure SMS_41
表示分析中心站坐标观测方程的设计矩阵；tr(·)为迹运算函数；
逐一计算去除每个测站后的方差改变量
Figure SMS_42
，并对/>
Figure SMS_43
进行从大到小排序处理，满足如下关系式时剔除站点观测量：/>
Figure SMS_44
其中，ndel表示待剔除的异常站坐标数量；
同时，若下式成立：
Figure SMS_45
则将第idel个测站标记为异常站点，在站坐标综合中剔除。
在以上技术方案的基础上，优选的，所述步骤S4中，所述利用方差分量估计方法进行不同分析中心产品的权重分配具体包括：
利用等价权函数计算每个测站的权重，计算公式如下：
Figure SMS_46
其中，
Figure SMS_47
表示分析中心s中测站isit的权重，c0为常数，/>
Figure SMS_48
表示测站isit的坐标残差距离，利用每个测站的权重构成对角阵，并与分析中心的权阵相乘作为下一次计算的分析中心权阵；
计算剔除各个分析中心站坐标产品后的粗差，重新计算后验方差，记为
Figure SMS_49
，按下式更新观测数据的权重缩放因子：/>
Figure SMS_50
其中，
Figure SMS_51
为权重缩放因子。
本发明第二方面，公开一种多分析中心精密站坐标产品综合系统，所述系统包括：
约束消除模块：用于对多个分析中心提供的GNSS站坐标产品文件进行预处理，利用分析中心提供的先验信息对预处理得到的原始法方程系统进行初步先验约束消除，利用相似变换矩阵消除残余约束；
框架对齐模块：用于建立相似变换方程，采用非潮汐海洋、大气和水文负载产品对目标框架进行非线性改正，将消除先验约束之后的分析中心站坐标产品对齐到目标框架，并计算目标框架下的分析中心站坐标解；
产品综合模块：用于对于所有的分析中心站坐标产品，根据分析中心站坐标解建立站坐标综合模型，得到总的综合法方程系统，利用最小二乘原理得到综合解估计值和相似变换参数以及站坐标残差，实现多分析中心站坐标产品综合；
异常处理模块：用于利用站点间的相关性信息识别分析中心站坐标产品中的异常站点；
权重分配模块：用于利用方差分量估计方法进行各个分析中心产品的权重分配；
循环迭代模块：用于重复产品综合模块、异常处理模块和权重分配模块，直至各个分析中心的权重趋于稳定或达到预设的最大迭代次数，得到权重分配结果；
求解输出模块：用于权重分配结果求解总的综合法方程系统，输出综合站坐标。
本发明第三方面，公开一种电子设备，包括：至少一个处理器、至少一个存储器、通信接口和总线；
其中，所述处理器、存储器、通信接口通过所述总线完成相互间的通信；
所述存储器存储有可被所述处理器执行的程序指令，所述处理器调用所述程序指令，以实现如本发明第一方面所述的方法。
本发明第四方面，公开一种计算机可读存储介质，所述计算机可读存储介质存储计算机指令，所述计算机指令使计算机实现如本发明第一方面所述的方法。
本发明相对于现有技术具有以下有益效果：
1）本发明综合利用了基于站坐相关性的粗差探测、先验残余约束消除方法和框架对齐策略，获得高可靠性和基准连续的综合站坐标产品；并在框架对齐策略中采用非潮汐海洋、大气和水文负载产品对目标框架进行非线性改正，可以改善分析中心站坐标网形和目标框架的相似性，提高框架对齐精度。
2）本发明基于一个相似变换矩阵，利用数学变换对分析中心站坐标产品恢复的法方程进行处理，可以消除残余的先验约束对站坐标网形的扭曲，增强产品综合结果的稳定性。
3）本发明利用站点间方差-协方差信息识别分析中心站坐标产品中的异常站点，提高观测数据质量，并利用后验方差迭代更新不同分析中心产品的权重，提高站坐标综合结果的可靠性。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1为本发明的多分析中心精密站坐标产品综合方法流程图。
具体实施方式
下面将结合本发明实施方式，对本发明实施方式中的技术方案进行清楚、完整地描述，显然，所描述的实施方式仅仅是本发明一部分实施方式，而不是全部的实施方式。基于本发明中的实施方式，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施方式，都属于本发明保护的范围。
本发明的实施例中，均假设多个分析中心提供的站坐标产品包含有公共的测站，只有这样才能有效地建立观测方程，实现站坐标产品的综合。在实际中该条件也非常容易满足，因为从iGMAS/IGS分析中心来看，有超过10家分析中心向全球用户提供不少于100个全球跟踪站的坐标结果，且所有分析中心站坐标产品中都包含数十个核心站的坐标解。站坐标综合中的观测方程是建立在多个站点坐标的基础上，缺失的个别站点数据不影响站坐标综合中的平差处理。
请参阅图1，本发明提出一种多分析中心精密站坐标产品综合方法，所述方法包括：
S1、对多个分析中心提供的GNSS站坐标产品文件进行预处理，利用分析中心提供的先验信息对预处理得到的原始法方程系统进行初步先验约束消除，利用相似变换矩阵消除残余约束。
步骤S1具体包括如下分步骤：
S11、获取多个分析中心提供的GNSS站坐标产品文件并进行预处理。
具体的，获取分析中心站坐标文件和ITRF框架解文件，初步判定各精密站坐标产品的格式，读取分析中心站坐标数据，并统计精密站坐标产品中测站数、核心站，恢复参数估计的原始法方程系统。
S12、初步消除先验约束。
设添加先验约束的站坐标估计系统如下：
Figure SMS_52
（1）
公式（1）中，
Figure SMS_53
为站坐标先验值，/>
Figure SMS_54
为站坐标估计值，/>
Figure SMS_55
为先验约束矩阵，/>
Figure SMS_56
为后验方差-协方差阵，/>
Figure SMS_57
为未添加任何约束的信息矩阵；u是原始法方程系统对应的右向量；
利用分析中心提供的先验约束信息对原始法方程系统进行初步处理，部分消除先验约束，采用如下公式初步消除先验约束：
Figure SMS_58
（2）
公式（2）中，
Figure SMS_59
为没有完全消除先验约束的信息矩阵。
S13、消除残余约束。
考虑到分析中心在站坐标解算中同时估计了轨道、钟差、对流层、模糊度等参数，而在输出的站坐标法方程系统时将这些参数都预先进行了消除，因此提供的是一个约化法方程系统，而这些消除的参数引入的先验约束会被传递到约化法方程中，因此利用公式（2）并不能得到完全无约束的法方程系统，需要进一步进行残余约束的消除。
本发明采用下式进行残余约束的消除：
Figure SMS_60
（3）
公式（3）中，
Figure SMS_61
为单位矩阵，/>
Figure SMS_62
为无约束法方程系统的右矩阵，/>
Figure SMS_63
为无约束对应的坐标估计值；N为原始法方程系统中的信息矩阵，u为原始法方程系统对应的右向量；/>
Figure SMS_64
表示有n个数值1组成的列向量，n为跟踪站数，/>
Figure SMS_65
为克罗内克积，M为相似变换矩阵，M的表达式为：/>
Figure SMS_66
（4）
公式（4）中，x i 、y i 、z i 分别为测站i的X、Y和Z坐标分量。
通过公式（4），可将基准相关的所有信息（原点、尺度、方向等）都进行去除，方便进行多分析中心站点综合。
本发明基于一个相似变换矩阵M，利用数学变换对分析中心站坐标产品恢复的原始法方程系统进行处理，可以消除残余的先验约束对站坐标网形的扭曲，增强产品综合结果的稳定性。
S2、建立相似变换方程，采用非潮汐海洋、大气和水文负载产品对目标框架进行非线性改正，将消除先验约束之后的分析中心站坐标产品对齐到目标框架，并计算目标框架下的分析中心站坐标解。
所述步骤S2具体包括如下分步骤：
S21、计算站点在当前历元t的非潮汐海洋、大气和水文负载形变，将非潮汐海洋、大气和水文负载形变之和记为
Figure SMS_67
。
S22、以国际地球参考框架解为目标框架，在目标框架中引入非潮汐海洋、大气和水文负载，计算目标框架在当前历元的坐标。
目标框架提供的是参考历元的位置和速度，是一种线性框架，本发明在该线性框架的基础上，引入非潮汐的海洋、大气和水文负载，形成非线性的目标框架，以改善目标框架和分析中心站坐标网形的相似性。
该非线性的目标框架在当前历元（分析中心站坐标产品的参考历元）的坐标计算公式为：
Figure SMS_68
（5）
公式（5）中，t和t 0分别表示当前历元和参考历元；
Figure SMS_69
表示目标框架解在当前历元的位置；/>
Figure SMS_70
表示目标框架解在参考历元的位置；/>
Figure SMS_71
表示目标框架解在参考历元的速度。
本发明采用非潮汐海洋、大气和水文负载产品对目标框架进行非线性改正，可以改善分析中心站坐标网形和目标框架的相似性，从而提高框架对齐精度。
S23、计算分析中心站坐标产品中核心站点目标框架解的坐标值
Figure SMS_72
。
S24、对分析中心无约束的法方程系统添加统一约束：
Figure SMS_73
（6）
公式（6）中，
Figure SMS_74
，是由核心站坐标构成的约束矩阵，/>
Figure SMS_75
为核心站坐标构成的相似变换矩阵，/>
Figure SMS_76
为相应的转置矩阵，/>
Figure SMS_77
为目标框架下的分析中心站坐标解，s表示分析中心编号。
S25、计算目标框架下的分析中心站坐标解
Figure SMS_78
。
基于公式（3）和（6），计算得到目标框架下的分析中心站坐标解，记为
Figure SMS_79
。
S3、对于所有的分析中心站坐标产品，根据分析中心站坐标解建立站坐标综合模型，得到总的综合法方程系统，利用最小二乘原理得到综合解估计值和相似变换参数以及站坐标残差。
根据分析中心站坐标解
Figure SMS_80
，建立如下站坐标综合模型：
Figure SMS_81
（7）
公式（7）中，s表示分心中心编号，
Figure SMS_82
表示分析中心站坐标残差，/>
Figure SMS_83
表示综合站坐标解；/>
Figure SMS_84
表示相似变换中的平移参数向量；/>
Figure SMS_85
表示相似变换中的旋转参数矩阵；
Figure SMS_86
相似变换中的尺度参数标量。
本发明对所有分析中心站坐标产品建立公式如公式（7）的观测方程，得到总的综合法方程系统，利用最小二乘原理得到综合解估计值和相似变换参数以及站坐标残差，初步实现多分析中心站坐标产品综合。
S4、利用站点间的相关性信息识别分析中心站坐标产品中的异常站点，并利用方差分量估计方法进行不同分析中心产品的权重分配。
步骤S4具体包括如下分步骤：
S41、利用站点间方差-协方差信息识别分析中心站坐标产品中的异常站点。
步骤S41具体可划分为如下步骤：
S411、将每个分析中心的站坐标产品划分为一类观测数据，采用下式计算每个分析中心的后验方差：
Figure SMS_87
（8）
公式（8）中，
Figure SMS_88
表示后验方差；nsit表示分析中心s的站坐标产品中测站总数；设/>
Figure SMS_89
表示分析中心s中的跟踪站isit的坐标后验残差，/>
Figure SMS_90
表示分析中心s中的跟踪站k的坐标后验残差，/>
Figure SMS_91
表示站点isit和k之间的相关信息；/>
Figure SMS_92
表示综合法方程系统的信息矩阵逆；/>
Figure SMS_93
表示分析中心站坐标观测量的权阵；/>
Figure SMS_94
表示分析中心站坐标观测方程的设计矩阵；tr(·)为迹运算函数。
S412、逐一计算去除每个测站后的方差改变量，记去除第idel个测站后的方差改变量为
Figure SMS_95
，对/>
Figure SMS_96
进行从大到小的排序处理，满足如下关系式时剔除站点的观测量：/>
Figure SMS_97
（9）
公式（9）中，ndel表示待剔除的异常站坐标数量。
站点观测量之间的相关性使得异常观测的影响将传递给其它观测量，甚至出现为后验方差
Figure SMS_98
为负值的情况，针对这种情况，本发明逐一计算去除每个测站后的方差改变量，从而通过粗差探测确定异常站点。
S413、在步骤S412的基础上，若下式成立：
Figure SMS_99
（10）
则将第idel个测站标记为异常站点，从站坐标综合中剔除。
S42、利用等价权函数计算每个测站的权重。
设
Figure SMS_100
表示分析中心s中测站isit的权重，权重计算公式如下：
Figure SMS_101
（11）
公式（11）中，
Figure SMS_102
表示分析中心s中测站isit的坐标残差距离，c 0 为经验常数，利用每个测站的权重构成对角阵，并与分析中心的权阵相乘作为下一次计算的分析中心权阵；
计算剔除各个分析中心站坐标产品后的粗差，重新计算后验方差，记为
Figure SMS_103
，按下式更新观测数据的权重缩放因子：/>
Figure SMS_104
（12）
其中，
Figure SMS_105
为权重缩放因子。
S5、重复以上步骤S3~S4，直至分析中心的权重趋于稳定或达到预设的最大迭代次数。
具体的，当前后两次迭代的分析中心权重缩放因子
Figure SMS_106
小于预设阈值0.1，则表示分析中心的权重趋于稳定，结束迭代；或者当迭代次数达到预设的最大迭代次数5次，结束迭代，得到权重分配结果。
S6、基于权重分配结果求解总的综合法方程系统，输出综合站坐标。
基于迭代得到的各个分析中心站坐标产品的权重，计算叠加之后的法方程系统，得到综合的站坐标估计值。
本发明在站坐标产品综合中综合利用了先验残余约束消除策略、框架对齐策略和基于站坐标相关性的粗差探测，可以获得高可靠性和基准连续的综合站坐标产品。其中，先验残余约束消除方法通过利用数学变换对分析中心站坐标产品恢复的法方程进行处理，消除残余的先验约束对站坐标网形的扭曲；框架对齐策略提高在目标框架中引入非线性站点形变改正，来改善分析中心站坐标网形和目标框架的相似性，提高框架对齐精度；基于站坐相关性的粗差探测利用站点间方差-协方差信息识别分析中心站坐标产品中的异常站点，并利用方差分量估计方法确定不同分析中心产品的权重，最终解决多分析中心精密站坐标产品综合中的技术难点。
与上述方法实施例相对应，本发明还提出一种多分析中心精密站坐标产品综合系统，所述系统包括：
约束消除模块：用于对多个分析中心提供的GNSS站坐标产品文件进行预处理，利用分析中心提供的先验信息对预处理得到的原始法方程系统进行初步先验约束消除，利用相似变换矩阵消除残余约束；
框架对齐模块：用于采用非潮汐海洋、大气和水文负载产品对目标框架进行非线性改正，将消除先验约束之后的分析中心站坐标产品对齐到目标框架，并计算目标框架下的分析中心站坐标解；
产品综合模块：用于对于所有的分析中心站坐标产品，根据分析中心站坐标解建立站坐标综合模型，得到总的综合法方程系统，利用最小二乘原理得到综合解估计值和相似变换参数以及站坐标残差，实现多分析中心站坐标产品综合；
异常处理模块：用于利用站点间的相关性信息识别分析中心站坐标产品中的异常站点；
权重分配模块：用于利用方差分量估计方法进行各个分析中心产品的权重分配；
循环迭代模块：用于重复产品综合模块、异常处理模块和权重分配模块，直至各个分析中心的权重趋于稳定或达到预设的最大迭代次数，得到权重分配结果；
求解输出模块：用于基于权重分配结果求解总的综合法方程系统，输出综合站坐标。
以上系统实施例和方法实施例是一一对应的，系统实施例简述之处请参阅方法实施例即可。
本发明还公开一种电子设备，包括：至少一个处理器、至少一个存储器、通信接口和总线；其中，所述处理器、存储器、通信接口通过所述总线完成相互间的通信；所述存储器存储有可被所述处理器执行的程序指令，所述处理器调用所述程序指令，以实现本发明前述的方法。
本发明还公开一种计算机可读存储介质，所述计算机可读存储介质存储计算机指令，所述计算机指令使所述计算机实现本发明实施例所述方法的全部或部分步骤。所述存储介质包括：U盘、移动硬盘、只读存储器ROM、随机存取存储器RAM、磁碟或者光盘等各种可以存储程序代码的介质。
以上所描述的系统实施例仅仅是示意性的，其中所述作为分离部件说明的单元可以是或者也可以不是物理上分开的，作为单元显示的部件可以是或者也可以不是物理单元，即可以分布到多个网络单元上。本领域普通技术人员在不付出创造性的劳动的情况下，可以根据实际的需要选择其中的部分或者全部模块来实现本实施例方案的目的。
以上所述仅为本发明的较佳实施方式而已，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (8)

1.一种多分析中心精密站坐标产品综合方法，其特征在于，所述方法包括：
S1、对多个分析中心提供的GNSS站坐标产品文件进行预处理，利用分析中心提供的先验信息对预处理得到的原始法方程系统进行初步先验约束消除，利用相似变换矩阵消除残余约束；
S2、采用非潮汐海洋、大气和水文负载产品对目标框架进行非线性改正，将消除先验约束之后的分析中心站坐标产品对齐到目标框架，并计算目标框架下的分析中心站坐标解；所述步骤S2具体包括如下分步骤：
计算站点在当前历元t的非潮汐海洋、大气和水文负载形变，非潮汐海洋、大气和水文负载形变记为
Figure QLYQS_1
；
在目标框架中引入非潮汐海洋、大气和水文负载，计算目标框架在当前历元的坐标，计算公式为：
Figure QLYQS_2
；
t和t 0 分别表示当前历元和参考历元；
Figure QLYQS_3
表示目标框架解在当前历元的位置；/>
Figure QLYQS_4
表示目标框架解在参考历元的位置；/>
Figure QLYQS_5
表示目标框架解在参考历元的速度；
计算分析中心站坐标产品中核心站点目标框架解的坐标值
Figure QLYQS_6
；
对分析中心无约束的法方程系统添加统一约束：
Figure QLYQS_7
；
其中，
Figure QLYQS_8
是由核心站坐标构成的约束矩阵，/>
Figure QLYQS_9
为核心站坐标构成的相似变换矩阵，/>
Figure QLYQS_10
为相应的转置矩阵，/>
Figure QLYQS_11
为目标框架下的分析中心站坐标解，s表示分析中心编号；
计算目标框架下的分析中心站坐标解
Figure QLYQS_12
；
S3、对于所有的分析中心站坐标产品，根据分析中心站坐标解建立站坐标综合模型，得到总的综合法方程系统，利用最小二乘原理得到综合解估计值和相似变换参数以及站坐标残差，实现多分析中心站坐标产品综合；
S4、利用站点间的相关性信息识别分析中心站坐标产品中的异常站点，并利用方差分量估计方法进行各个分析中心产品的权重分配；
所述利用站点间的相关性信息识别分析中心站坐标产品中的异常站点具体包括：
将每个分析中心的站坐标产品划分为一类观测数据，采用下式计算每个分析中心的后验方差：
Figure QLYQS_13
；
其中，
Figure QLYQS_14
表示后验方差；nsit表示分析中心s的站坐标产品中的站点总数；设
Figure QLYQS_15
表示分析中心s中的跟踪站isit的坐标后验残差，/>
Figure QLYQS_16
表示分析中心s中的跟踪站k的坐标后验残差，/>
Figure QLYQS_17
表示站点isit和k之间的相关信息；/>
Figure QLYQS_18
表示综合法方程系统的信息矩阵逆；/>
Figure QLYQS_19
表示分析中心站坐标观测量的权阵；/>
Figure QLYQS_20
表示分析中心站坐标观测方程的设计矩阵；tr(·)为迹运算函数；
逐一计算去除每个站点后的方差改变量
Figure QLYQS_21
，并对/>
Figure QLYQS_22
进行从大到小排序处理，满足如下关系式时剔除站点观测量：
Figure QLYQS_23
；
其中，ndel表示待剔除的异常站坐标数量；
同时，若下式成立：
Figure QLYQS_24
；
则将第idel个站点标记为异常站点，在站坐标综合中剔除；
S5、重复以上步骤S3~S4，直至各个分析中心的权重趋于稳定或达到预设的最大迭代次数，得到权重分配结果；
S6、基于权重分配结果求解总的综合法方程系统，输出综合站坐标。
2.根据权利要求1所述的多分析中心精密站坐标产品综合方法，其特征在于，所述预处理具体包括：初步判定各精密站坐标产品的格式，并统计精密站坐标产品中站点数、核心站数以及存储原始法方程系统信息。
3.根据权利要求2所述的多分析中心精密站坐标产品综合方法，其特征在于，所述利用相似变换矩阵消除残余约束具体包括：
设添加先验约束的站坐标估计系统如下：
Figure QLYQS_25
；
其中，
Figure QLYQS_26
为站坐标先验值，/>
Figure QLYQS_27
为站坐标估计值，/>
Figure QLYQS_28
为先验约束矩阵，/>
Figure QLYQS_29
为后验方差-协方差阵，/>
Figure QLYQS_30
为未添加任何约束的信息矩阵，/>
Figure QLYQS_31
为原始法方程系统对应的右向量；
初步先验约束消除采用的公式如下：
Figure QLYQS_32
；
其中，
Figure QLYQS_33
为没有完全消除先验约束的信息矩阵；
采用下式进行残余约束的消除：
Figure QLYQS_34
；
其中，
Figure QLYQS_35
为单位矩阵，/>
Figure QLYQS_36
为无约束法方程系统的右矩阵，/>
Figure QLYQS_37
为无约束对应的坐标估计值；N为原始法方程系统中的信息矩阵，/>
Figure QLYQS_38
表示有n个数值1组成的列向量，n为跟踪站数，
Figure QLYQS_39
为克罗内克积，M为相似变换矩阵，M的表达式为：
Figure QLYQS_40
；
其中x i 、y i 、z i 分别为站点i的X、Y和Z坐标分量。
4.根据权利要求1所述的多分析中心精密站坐标产品综合方法，其特征在于，所述站坐标综合模型的公式为：
Figure QLYQS_41
；
其中，s表示分心中心编号，
Figure QLYQS_42
表示分析中心站坐标残差，/>
Figure QLYQS_43
表示综合站坐标解；
Figure QLYQS_44
表示相似变换中的平移参数向量；/>
Figure QLYQS_45
表示相似变换中的旋转参数矩阵；/>
Figure QLYQS_46
相似变换中的尺度参数标量。
5.根据权利要求1所述的多分析中心精密站坐标产品综合方法，其特征在于，所述步骤S4中，所述利用方差分量估计方法进行不同分析中心产品的权重分配具体包括：
利用等价权函数计算每个站点的权重，计算公式如下：
Figure QLYQS_47
；
其中，
Figure QLYQS_48
表示分析中心s中站点isit的权重，c0为常数，/>
Figure QLYQS_49
表示站点isit的坐标残差距离，利用每个站点的权重构成对角阵，并与分析中心的权阵相乘作为下一次计算的分析中心权阵；
计算剔除各个分析中心站坐标产品后的粗差，重新计算后验方差，记为
Figure QLYQS_50
，按下式更新观测数据的权重缩放因子：
Figure QLYQS_51
；
其中，
Figure QLYQS_52
为权重缩放因子。
6.使用权利要求1~5任一项所述方法的一种多分析中心精密站坐标产品综合系统，其特征在于，所述系统包括：
约束消除模块：用于对多个分析中心提供的GNSS站坐标产品文件进行预处理，利用分析中心提供的先验信息对预处理得到的原始法方程系统进行初步先验约束消除，利用相似变换矩阵消除残余约束；
框架对齐模块：用于建立相似变换方程，采用非潮汐海洋、大气和水文负载产品对目标框架进行非线性改正，将消除先验约束之后的分析中心站坐标产品对齐到目标框架，并计算目标框架下的分析中心站坐标解；
产品综合模块：用于对于所有的分析中心站坐标产品，根据分析中心站坐标解建立站坐标综合模型，得到总的综合法方程系统，利用最小二乘原理得到综合解估计值和相似变换参数以及站坐标残差，实现多分析中心站坐标产品综合；
异常处理模块：用于利用站点间的相关性信息识别分析中心站坐标产品中的异常站点；
权重分配模块：用于利用方差分量估计方法进行各个分析中心产品的权重分配；
循环迭代模块：用于重复产品综合模块、异常处理模块和权重分配模块，直至各个分析中心的权重趋于稳定或达到预设的最大迭代次数，得到权重分配结果；
求解输出模块：用于基于权重分配结果求解总的综合法方程系统，输出综合站坐标。
7.一种电子设备，其特征在于，包括：至少一个处理器、至少一个存储器、通信接口和总线；
其中，所述处理器、存储器、通信接口通过所述总线完成相互间的通信；
所述存储器存储有可被所述处理器执行的程序指令，所述处理器调用所述程序指令，以实现如权利要求1~5任一项所述的方法。
8.一种计算机可读存储介质，其特征在于，所述计算机可读存储介质存储计算机指令，所述计算机指令使计算机实现如权利要求1~5任一项所述的方法。
CN202310480488.1A 2023-04-28 2023-04-28 一种多分析中心精密站坐标产品综合方法及系统 Active CN116204756B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202310480488.1A CN116204756B (zh) 2023-04-28 2023-04-28 一种多分析中心精密站坐标产品综合方法及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202310480488.1A CN116204756B (zh) 2023-04-28 2023-04-28 一种多分析中心精密站坐标产品综合方法及系统

Publications (2)

Publication Number Publication Date
CN116204756A CN116204756A (zh) 2023-06-02
CN116204756B true CN116204756B (zh) 2023-07-07

Family

ID=86508011

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202310480488.1A Active CN116204756B (zh) 2023-04-28 2023-04-28 一种多分析中心精密站坐标产品综合方法及系统

Country Status (1)

Country Link
CN (1) CN116204756B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116859421B (zh) * 2023-06-27 2024-03-15 国汽大有时空科技(安庆)有限公司 一种多参考框架的定位服务方法及装置
CN116626724B (zh) * 2023-07-24 2023-10-10 齐鲁空天信息研究院 一种基于数字广播的ssr信息传输与评估方法
CN117388872B (zh) * 2023-09-05 2024-03-19 武汉大学 一种北斗地基增强系统参考站坐标框架维持方法和系统

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2012021760A2 (en) * 2010-08-12 2012-02-16 The Government Of The United States Of America As Represented By The Secretary Of The Navy Improved orbit covariance estimation and analysis (ocean) system and method
CN110081909A (zh) * 2019-05-22 2019-08-02 北京中交华安科技有限公司 基于全球定位控制点坐标的车载移动测量系统检校方法
WO2019205299A1 (zh) * 2018-04-27 2019-10-31 中国农业大学 视觉测量系统结构参数标定和仿射坐标系构建方法与系统
CN110398753A (zh) * 2019-06-28 2019-11-01 武汉大学 Gnss测站坐标时间序列周期性探测方法及系统

Family Cites Families (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8760343B2 (en) * 2009-11-17 2014-06-24 Topcon Positioning Systems, Inc. Detection and correction of anomalous measurements and ambiguity resolution in a global navigation satellite system receiver
CN103591891B (zh) * 2013-11-20 2015-04-29 天津大学 室内空间测量定位系统的精密控制场精度溯源方法
CN105572703B (zh) * 2015-12-17 2016-09-28 武汉大学 一种gps时间序列广义共模误差提取方法
CN108415039B (zh) * 2018-01-25 2022-03-15 中国人民解放军63921部队 iGMAS多分析中心多卫星系统精密轨道产品综合方法
CN108536648B (zh) * 2018-03-30 2021-07-06 武汉大学 基于多超声波传感器的局部放电非线性模型转换求解与优化方法
CN108845340A (zh) * 2018-06-01 2018-11-20 浙江亚特电器有限公司 基于gnss-rtk的定位方法
US11138465B2 (en) * 2019-12-10 2021-10-05 Toyota Research Institute, Inc. Systems and methods for transforming coordinates between distorted and undistorted coordinate systems
FR3106241B1 (fr) * 2020-01-10 2023-04-28 Alstom Transp Tech Système et procédé de gestion de l’énergie d’alimentation d’un véhicule de transport, et véhicule de transport correspondant
CN112799101A (zh) * 2021-01-29 2021-05-14 华东师范大学 一种构建gnss区域大地参考框架的方法
CN113848577A (zh) * 2021-08-19 2021-12-28 中国能源建设集团江苏省电力设计院有限公司 一种基于动态分区的大规模gnss网并行解算方法及系统
CN115982564A (zh) * 2022-12-05 2023-04-18 广东电网有限责任公司 一种区域电离层电子密度计算方法、装置及计算机设备

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2012021760A2 (en) * 2010-08-12 2012-02-16 The Government Of The United States Of America As Represented By The Secretary Of The Navy Improved orbit covariance estimation and analysis (ocean) system and method
WO2019205299A1 (zh) * 2018-04-27 2019-10-31 中国农业大学 视觉测量系统结构参数标定和仿射坐标系构建方法与系统
CN110081909A (zh) * 2019-05-22 2019-08-02 北京中交华安科技有限公司 基于全球定位控制点坐标的车载移动测量系统检校方法
CN110398753A (zh) * 2019-06-28 2019-11-01 武汉大学 Gnss测站坐标时间序列周期性探测方法及系统

Also Published As

Publication number Publication date
CN116204756A (zh) 2023-06-02

Similar Documents

Publication Publication Date Title
CN116204756B (zh) 一种多分析中心精密站坐标产品综合方法及系统
Banerjee et al. Nearest neighbour distributions: New statistical measures for cosmological clustering
CN107765275B (zh) 广域差分定位方法、装置、终端及计算机可读存储介质
CN110568459B (zh) 基于igs和cors站的区域电离层tec实时监测方法
Brockmann et al. An improved model of the Earth’s static gravity field solely derived from reprocessed GOCE data
Liu et al. The cooperative IGS RT-GIMs: A reliable estimation of the global ionospheric electron content distribution in real time
Maksyutov et al. A high-resolution inverse modelling technique for estimating surface CO 2 fluxes based on the NIES-TM–FLEXPART coupled transport model and its adjoint
Mohammed et al. An assessment of static precise point positioning using GPS only, GLONASS only, and GPS plus GLONASS
CN112902989B (zh) 数据误差标定方法、装置、电子设备和存储介质
Rebischung et al. Recent results from the IGS terrestrial frame combinations
CN107422342A (zh) Gnss卫星钟差实时估计质量控制方法
Raymund et al. Model‐assisted ionospheric tomography: A new algorithm
CN115390095A (zh) 一种获取电离层延迟的方法、装置及介质
Li et al. Estimation and analysis of BDS2 and BDS3 differential code biases and global ionospheric maps using BDS observations
CN114662059A (zh) 一种海上卫星大地坐标的高程拟合方法与装置
Liu et al. Cloud-based single-frequency snapshot RTK positioning
Mirmohammadian et al. Multi-GNSS-Weighted interpolated tropospheric delay to improve long-baseline RTK positioning
Banville et al. Wide-area grid-based slant ionospheric delay corrections for precise point positioning
CN106569232A (zh) 空间信号精度评估方法及系统
Hobiger et al. Observation level combination of SLR and VLBI with c5++: a case study for TIGO
Favenza et al. A machine learning approach to GNSS scintillation detection: Automatic soft inspection of the events
Zhou et al. A robust trend estimator for GNSS time series in the presence of complex periodicity and its evaluation on multi-source products of IGS and IGMAS
Zhang et al. Centimeter-level positioning by instantaneous lidar-aided GNSS ambiguity resolution
Hunegnaw et al. Status of TIGA activities at the British Isles continuous GNSS Facility and the University of Luxembourg
Gourine Use of Starlette and LAGEOS-1&-2 laser measurements for determination and analysis of stations coordinates and EOP time series

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