CN114662268A - An Improved GNSS Network Sequential Adjustment Calculation Method - Google Patents

An Improved GNSS Network Sequential Adjustment Calculation Method Download PDF

Info

Publication number
CN114662268A
CN114662268A CN202111290536.8A CN202111290536A CN114662268A CN 114662268 A CN114662268 A CN 114662268A CN 202111290536 A CN202111290536 A CN 202111290536A CN 114662268 A CN114662268 A CN 114662268A
Authority
CN
China
Prior art keywords
adjustment
parameter
value
fuzzy
error
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.)
Granted
Application number
CN202111290536.8A
Other languages
Chinese (zh)
Other versions
CN114662268B (en
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.)
Guangzhou Urban Planning Survey and Design Institute
Original Assignee
Guangzhou Urban Planning Survey and Design Institute
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 Guangzhou Urban Planning Survey and Design Institute filed Critical Guangzhou Urban Planning Survey and Design Institute
Priority to CN202111290536.8A priority Critical patent/CN114662268B/en
Publication of CN114662268A publication Critical patent/CN114662268A/en
Application granted granted Critical
Publication of CN114662268B publication Critical patent/CN114662268B/en
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
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • 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
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/04—Constraint-based CAD

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Mathematical Optimization (AREA)
  • Data Mining & Analysis (AREA)
  • General Engineering & Computer Science (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Pure & Applied Mathematics (AREA)
  • Computer Hardware Design (AREA)
  • Operations Research (AREA)
  • Evolutionary Computation (AREA)
  • Algebra (AREA)
  • Geometry (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Data Exchanges In Wide-Area Networks (AREA)
  • Position Fixing By Use Of Radio Waves (AREA)

Abstract

The invention discloses an improved GNSS network sequential adjustment calculation method, which comprises the steps of establishing a first group of error equations and a second group of error equations of a sequential adjustment model, and performing single adjustment on the first group of error equations to obtain a parameter correction value matrix and a parameter covariance matrix; calculating according to the parameter correction value matrix to obtain a first adjustment value of the unknown parameter; calculating a diagonal value according to the parameter covariance matrix to obtain a medium error of an unknown parameter; determining a fuzzy center and a fuzzy amplitude of a public parameter correction value according to a fuzzy theory, and solving to obtain a parameter second correction value according to a constructed constrained adjustment function model; calculating a second adjustment value of the unknown parameter according to the second parameter correction value and the first adjustment value; and calculating to obtain the coordinates of each GNSS network point according to the second adjustment value. The parameter estimation distortion caused by gross error can be effectively weakened, error accumulation is reduced, and the resolving precision is improved.

Description

一种改进的GNSS网序贯平差计算方法An Improved GNSS Network Sequential Adjustment Calculation Method

技术领域technical field

本发明涉及卫星定位技术领域，尤其涉及一种改进的GNSS网序贯平差计算 方法、装置、存储介质及终端设备。The present invention relates to the technical field of satellite positioning, and in particular, to an improved GNSS network sequential adjustment calculation method, device, storage medium and terminal device.

背景技术Background technique

随着全球导航卫星系统(GNSS)飞速发展，利用GNSS建立各等级控制网 在各行各业得到了广泛应用，多个国家、地区、行业纷纷建立了满足自身要求的 GNSS基准站网。利用GNSS建立各等级基准站网均采用相对定位技术，即确定 测量点间的相对位置关系，将这种点间的相对位置量称为基线向量坐标，由基线 向量组成的观测网称为基线向量网，GNSS网平差就是以GNSS基线向量为观测 值进行平差计算获得各GNSS网点的坐标并进行精度评定的过程。With the rapid development of the Global Navigation Satellite System (GNSS), the use of GNSS to establish control networks at various levels has been widely used in all walks of life. Many countries, regions and industries have established GNSS reference station networks that meet their own requirements. The relative positioning technology is used to establish the base station network of each level by using GNSS, that is, the relative position relationship between the measurement points is determined. GNSS network adjustment is the process of using the GNSS baseline vector as the observed value to carry out adjustment calculation to obtain the coordinates of each GNSS network point and to evaluate the accuracy.

在进行大规模GNSS网整体解算时，现有技术一般采用序贯平差估计完成 GNSS网点的坐标解算及其精度评定。将整个GNSS网分为若干个子网，将各个 子网独立解算，得到松约束下的参数估值及其协方差阵，然后将各子网联合处理。 序贯平差估计利用前期平差结果与当期观测样本获得与整体平差结果相同的最 优解。When performing the overall calculation of a large-scale GNSS network, the prior art generally adopts sequential adjustment estimation to complete the coordinate calculation and accuracy evaluation of the GNSS network points. The entire GNSS network is divided into several sub-networks, and each sub-network is solved independently to obtain the parameter estimation and its covariance matrix under the loose constraints, and then the sub-networks are jointly processed. Sequential adjustment estimation uses previous adjustment results and current observation samples to obtain the same optimal solution as the overall adjustment results.

在GNSS观测值中，由于受到观测信号、传播路径和接收机等部分的影响， 观测值中不可避免地带有粗差，但是由于现有技术是通过法方程叠加来完成各子 网间的联合处理，其本质是最小二乘，对粗差不具有抵抗力，当观测样本含有粗 差时，无法求出准确得平差值。In the GNSS observation value, due to the influence of the observation signal, propagation path and receiver, etc., the observation value inevitably contains gross errors, but because the existing technology is to complete the joint processing between the sub-networks through the superposition of normal equations. , whose essence is least squares, is not resistant to gross errors. When the observed sample contains gross errors, it is impossible to obtain an accurate adjustment value.

发明内容SUMMARY OF THE INVENTION

本发明实施例提供一种改进的GNSS网序贯平差计算方法，能够减少观测样 本粗差对于后续平差估值的影响，减小误差积累效应，输出准确的平差值。The embodiment of the present invention provides an improved GNSS network sequential adjustment calculation method, which can reduce the influence of gross errors of observed samples on subsequent adjustment estimates, reduce the effect of error accumulation, and output accurate adjustment values.

本发明实施例提供一种改进的GNSS网序贯平差计算方法，所述方法包括：The embodiment of the present invention provides an improved GNSS network sequential adjustment calculation method, and the method includes:

通过定义GNSS网中前期子网独立测站坐标参数，子网间公共测站坐标参数， 后期子网独立测站坐标参数，根据各测站间基线向量的几何关系建立前后期平差 模型的第一组误差方程和第二组误差方程；By defining the coordinate parameters of the independent stations of the early sub-network in the GNSS network, the coordinate parameters of the public stations between the sub-networks, and the coordinate parameters of the independent stations of the sub-network in the later period, the first step of the adjustment model before and after the adjustment model is established according to the geometric relationship of the baseline vectors between the stations. a set of error equations and a second set of error equations;

对所述第一组平差误差方程单独进行平差，得到参数改正值矩阵和参数协方 差矩阵；Adjusting the first group of adjustment error equations independently to obtain a parameter correction value matrix and a parameter covariance matrix;

根据所述参数改正值矩阵计算得到未知参数的第一次平差值；Calculate the first adjustment value of the unknown parameter according to the parameter correction value matrix;

对所述参数方差-协方差矩阵取对角值计算得到未知参数的中误差；Taking the diagonal value of the parameter variance-covariance matrix to calculate the median error of the unknown parameter;

将所述公共参数的第一次平差值作为第二组平差的近似值，代入所述第二组 平差误差方程中，计算新的常数项，定义新的观测信息改正数为V′2，得到新的误 差方程；Take the first adjustment value of the common parameter as the approximate value of the second group of adjustment, and substitute it into the second group of adjustment error equations, calculate a new constant term, and define the new observation information correction number as V' 2 , get a new error equation;

根据模糊理论，以所述第一次平差值和所述中误差确定公共参数改正值的模 糊中心和模糊幅度，根据所述常数项、所述模糊中心和所述模糊幅度构建平差函 数约束模型；According to fuzzy theory, the first adjustment value and the intermediate error are used to determine the fuzzy center and the fuzzy amplitude of the correction value of the common parameter, and the adjustment function constraint is constructed according to the constant term, the fuzzy center and the fuzzy amplitude Model;

对所述构建约束平差函数模型求解得到参数的第二改正值；The second correction value of the parameter is obtained by solving the construction constraint adjustment function model;

根据所述第二改正值和所述第一次平差值计算得到未知参数的第二次平差 值；Calculate the second adjustment value of the unknown parameter according to the second correction value and the first adjustment value;

根据所述第二次平差值计算获得各GNSS网点的坐标。The coordinates of each GNSS network point are obtained by calculating according to the second adjustment value.

优选地，所述第一组误差方程为V1＝A11Xa+A12Xb-f1；Preferably, the first set of error equations is V 1 =A 11 X a +A 12 X b -f 1 ;

所述第二组误差方程为V2＝B22Xb+B23Y-f2；The second set of error equations is V 2 =B 22 X b +B 23 Yf 2 ;

其中，V1为第一次观测值改正数，A11和A12为第一次平差模型系数矩阵，V2为第二次观测值改正数，B22和B23为第二次平差模型系数矩阵，Xa为所述前期子 网独立站点坐标参数，Xb为所述子网间公共站点坐标参数，Y为所述后期子网新 增站点坐标参数，f1为第一次平差误差方程的常数项，

Figure BDA0003334540950000031
f2为第二次平差误差方程的常数项，
Figure BDA0003334540950000032
L1和L2分别为第一 次平差观测值和第二次平差观测值，
Figure BDA0003334540950000033
和Y0为未知参数第一次平差解算时 所取的近似值。Among them, V 1 is the correction number of the first observation, A 11 and A 12 are the coefficient matrix of the first adjustment model, V 2 is the correction number of the second observation, and B 22 and B 23 are the second adjustment. Model coefficient matrix, X a is the coordinate parameter of the independent site of the sub-network in the early stage, X b is the coordinate parameter of the public site between the sub-networks, Y is the coordinate parameter of the newly added site of the sub-network in the later stage, and f 1 is the first flattening parameter. the constant term of the difference error equation,
Figure BDA0003334540950000031
f 2 is the constant term of the second adjustment error equation,
Figure BDA0003334540950000032
L 1 and L 2 are the first adjustment observations and the second adjustment observations, respectively,
Figure BDA0003334540950000033
and Y 0 are the approximations taken during the first adjustment of the unknown parameters.

作为一种优选方式，所述参数改正值矩阵为

Figure BDA0003334540950000034
As a preferred way, the parameter correction value matrix is
Figure BDA0003334540950000034

所述参数协方差矩阵为

Figure BDA0003334540950000035
The parameter covariance matrix is
Figure BDA0003334540950000035

其中，

Figure BDA0003334540950000036
和
Figure BDA0003334540950000037
为未知参数的改正数，P1为观测值权矩阵，
Figure BDA0003334540950000038
为单位 权方差,r为第一次平差时多余观测数。。in,
Figure BDA0003334540950000036
and
Figure BDA0003334540950000037
is the correction number of the unknown parameter, P 1 is the observation weight matrix,
Figure BDA0003334540950000038
is the unit weight variance, and r is the number of excess observations in the first adjustment. .

优选地，所述第一次平差值为

Figure BDA0003334540950000039
Preferably, the first adjustment value is
Figure BDA0003334540950000039

优选地，所述公共参数协方差为

Figure BDA00033345409500000310
为验后单位权方差。Preferably, the common parameter covariance is
Figure BDA00033345409500000310
is the posterior unit weight variance.

作为一种优选方式，所述将所述公共参数第一次平差值作为第二次平差时的 近似值，代入所述第二组平差误差方程中，计算新的常数项，定义新的观测信息 改正数为V′2，得到新的误差方程，具体包括：As a preferred way, the first adjustment value of the common parameter is used as an approximate value during the second adjustment, and is substituted into the second set of adjustment error equations, a new constant term is calculated, and a new The correction number of the observation information is V′ 2 , and a new error equation is obtained, which includes:

将所述第一次平差值

Figure BDA00033345409500000311
作为第二次平差时的近似值，代入所述第二组平差 误差方程中，计算得到新的常数项l2，定义新的观测信息改正数为V′2，得到新的 误差方程；the first adjustment value
Figure BDA00033345409500000311
As the approximate value of the second adjustment, it is substituted into the second set of adjustment error equations, and a new constant term l 2 is obtained by calculation, and the new observation information correction number is defined as V′ 2 , and a new error equation is obtained;

其中，

Figure BDA00033345409500000312
in,
Figure BDA00033345409500000312

优选地，所述根据模糊理论，以所述公共参数第一次平差值和所述中误差确 定参数改正值的模糊中心和模糊幅度，根据所述常数项、所述模糊中心和所述模 糊幅度构建平差函数约束模型，具体包括：Preferably, according to the fuzzy theory, the first adjustment value of the common parameter and the intermediate error are used to determine the fuzzy center and the fuzzy amplitude of the parameter correction value, and the fuzzy center and the fuzzy amplitude are determined according to the constant term, the fuzzy center and the fuzzy Amplitude builds an adjustment function constraint model, including:

以公共参数的第一次平差值

Figure BDA0003334540950000041
作为参数的模糊中心，则参数改正值的模糊中 心
Figure BDA0003334540950000042
为第二次平差时参数所取近似值。以所述中误差的3倍值为 模糊幅度Δ前；Take the first adjustment value of the common parameter
Figure BDA0003334540950000041
As the fuzzy center of the parameter, the fuzzy center of the parameter correction value
Figure BDA0003334540950000042
Approximate values for the parameters in the second adjustment. Taking 3 times of the middle error as the blur amplitude Δbefore ;

根据隶属函数、所述模糊中心和所述模糊幅度构建所述平差函数模型：The adjustment function model is constructed according to the membership function, the fuzzy center and the fuzzy magnitude:

Figure BDA0003334540950000043
Figure BDA0003334540950000043

其中，x″b和y′为第二次平差时参数的改正数，μA(x″b)为x″b的隶属函数，

Figure BDA0003334540950000044
Among them, x″ b and y′ are the correction numbers of the parameters during the second adjustment, μ A (x″ b ) is the membership function of x″ b ,
Figure BDA0003334540950000044

进一步地，所述对所述构建平差约束函数模型求解得到参数的第二改正值， 具体包括：Further, the second correction value of the parameter obtained by solving the construction adjustment constraint function model specifically includes:

对观测残差平方和取最小值的同时，x″b的隶属函数μA(x″b)取最大值，，得 到准则函数

Figure BDA0003334540950000045
While taking the minimum value of the sum of squares of the observed residuals, the membership function μ A (x″ b ) of x″ b takes the maximum value, and the criterion function is obtained
Figure BDA0003334540950000045

根据所述模糊幅度建立算子

Figure BDA0003334540950000046
Build an operator based on the blurring magnitude
Figure BDA0003334540950000046

根据所述算子将所述准则函数转化为准则函数矩阵

Figure BDA0003334540950000047
Transform the criterion function into a criterion function matrix according to the operator
Figure BDA0003334540950000047

对所述准则函数矩阵求偏导并令其等于0，计算得到参数的第二改正值

Figure BDA0003334540950000048
Calculate the partial derivative of the criterion function matrix and make it equal to 0, and calculate the second corrected value of the parameter
Figure BDA0003334540950000048

其中，τ为任一数值，0<τ<1，W＝diag[w1 w2 … wt]，Pi为观测值权阵，

Figure BDA0003334540950000049
为观测值残差，n＝1，2，3…，t＝1，2，3…，j＝1，2，...，t，Vxb＝x″b-xb前，表 示参数改正值与其先验模糊中心的偏差。Among them, τ is any value, 0<τ<1, W=diag[w 1 w 2 … w t ], P i is the observation weight matrix,
Figure BDA0003334540950000049
is the residual error of the observation value, n=1, 2, 3..., t=1, 2, 3..., j=1, 2,..., t, V xb = x″ b -x b before , it means the parameter correction The deviation of the value from its prior fuzzy center.

优选地，所述根据所述第二改正值和所述第一次平差值计算得到未知参数的 第二次平差值，具体为：Preferably, calculating the second adjustment value of the unknown parameter according to the second correction value and the first adjustment value, specifically:

将所述第二改正值

Figure BDA0003334540950000051
和所述第 一次平差值
Figure BDA0003334540950000052
代入平差值计算公式计算第二次平差值；the second correction value
Figure BDA0003334540950000051
and the first adjustment value
Figure BDA0003334540950000052
Substitute into the adjustment value calculation formula to calculate the second adjustment value;

所述平差值计算公式为

Figure BDA0003334540950000053
The calculation formula of the adjustment value is:
Figure BDA0003334540950000053

本发明提供一种改进的GNSS网序贯平差计算方法，通过建立前后期平差模 型的第一组误差方程和第二组误差方程；对所述第一组误差方程单独进行平差， 得到参数改正值矩阵和参数协方差矩阵；根据所述参数改正值矩阵计算得到未知 参数的第一次平差值；对参数协方差矩阵取对角值计算得到公共参数的中误差； 根据模糊理论，以所述第一次平差值和所述中误差确定参数的模糊中心和模糊幅 度，根据所述常数项、所述模糊中心和所述模糊幅度构建平差函数约束模型；对 所述构建约束平差函数模型求解得到参数的第二改正值；根据所述第二改正值和 所述第一次平差值计算得到未知参数的第二次平差值；根据所述第二次平差值计 算获得各GNSS网点的坐标。当后期观测信息含有粗差时，可以有效削弱粗差带 来的参数估值扭曲，减少误差积累，提高解算精度。The present invention provides an improved GNSS network sequential adjustment calculation method. By establishing the first group of error equations and the second group of error equations of the adjustment models before and after the adjustment, and adjusting the first group of error equations separately, the result is obtained: The parameter correction value matrix and the parameter covariance matrix; the first adjustment value of the unknown parameter is calculated according to the parameter correction value matrix; the diagonal value of the parameter covariance matrix is calculated to obtain the intermediate error of the public parameter; according to the fuzzy theory, Determine the fuzzy center and fuzzy amplitude of the parameter with the first adjustment value and the intermediate error, and construct an adjustment function constraint model according to the constant term, the fuzzy center and the fuzzy amplitude; for the construction constraint The adjustment function model is solved to obtain the second correction value of the parameter; the second adjustment value of the unknown parameter is calculated according to the second correction value and the first adjustment value; according to the second adjustment value Calculate the coordinates of each GNSS network point. When the later observation information contains gross errors, it can effectively weaken the distortion of parameter estimation caused by gross errors, reduce the accumulation of errors, and improve the calculation accuracy.

附图说明Description of drawings

图1是本发明实施例提供的一种改进的GNSS网序贯平差计算方法的流程示 意图；Fig. 1 is the schematic flow sheet of a kind of improved GNSS network sequential adjustment calculation method that the embodiment of the present invention provides;

图2所示是本发明实施例提供的GNSS网的网型图；2 is a network diagram of a GNSS network provided by an embodiment of the present invention;

图3是本发明实施例提供的序贯最小二乘及约束序贯算法的数据示意图。FIG. 3 is a schematic data diagram of sequential least squares and constrained sequential algorithms provided by an embodiment of the present invention.

具体实施方式Detailed ways

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、 完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的 实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前 提下所获得的所有其他实施例，都属于本发明保护的范围。The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative work fall within the protection scope of the present invention.

本发明实施例提供一种改进的GNSS网序贯平差计算方法，参见图1所示， 是本发明实施例提供的一种改进的GNSS网序贯平差计算方法的流程示意图，包 括步骤S1～S9：An embodiment of the present invention provides an improved GNSS network sequential adjustment calculation method, as shown in FIG. 1 , which is a schematic flowchart of an improved GNSS network sequential adjustment calculation method provided by an embodiment of the present invention, including step S1 ~S9:

S1，S1,

通过定义GNSS网中前期子网独立测站坐标参数，子网间公共测站坐标参数， 后期子网独立测站坐标参数，根据各测站间基线向量的几何关系建立前后期平差 模型的第一组误差方程和第二组误差方程；By defining the coordinate parameters of the independent stations of the early sub-network in the GNSS network, the coordinate parameters of the public stations between the sub-networks, and the coordinate parameters of the independent stations of the sub-network in the later period, the first step of the adjustment model before and after the adjustment model is established according to the geometric relationship of the baseline vectors between the stations. a set of error equations and a second set of error equations;

S2，对所述第一组平差误差方程进行单独平差，得到参数改正值矩阵和参数 协方差矩阵；S2, carry out a separate adjustment to the first group of adjustment error equations to obtain a parameter correction value matrix and a parameter covariance matrix;

S3，根据所述参数改正值矩阵计算得到未知参数的第一次平差值；S3, calculating the first adjustment value of the unknown parameter according to the parameter correction value matrix;

S4,对所述参数协方差矩阵取对角值计算得到公共参数的中误差；S4, taking the diagonal value of the parameter covariance matrix to calculate the median error of the public parameter;

S5，将所述第一次平差值作为第二次平差的近似值，代入所述第二组误差方 程中，计算新的常数项，定义新的观测信息改正数为V′2，得到新的误差方程；S5, taking the first adjustment value as the approximate value of the second adjustment, and substituting it into the second set of error equations, calculating a new constant term, defining a new observation information correction number as V′ 2 , and obtaining a new the error equation;

S6，根据模糊理论，以所述公共参数第一次平差值和所述中误差确定参数的 模糊中心和模糊幅度，根据所述常数项、所述模糊中心和所述模糊幅度构建平差 函数约束模型；S6, according to the fuzzy theory, determine the fuzzy center and the fuzzy amplitude of the parameter with the first adjustment value of the common parameter and the intermediate error, and construct an adjustment function according to the constant term, the fuzzy center and the fuzzy amplitude Constraint model;

S7，对所述构建平差约束函数模型求解得到参数的第二改正值；S7, the second correction value of the parameter is obtained by solving the construction adjustment constraint function model;

S8，根据所述第二改正值和所述第一次平差值计算得到未知参数的第二次平 差值；S8, calculates the second adjustment value of unknown parameter according to the second correction value and the first adjustment value;

S9，根据所述第二次平差值计算获得各GNSS网点的坐标。S9, calculate and obtain the coordinates of each GNSS network point according to the second adjustment value.

在本实施例具体实施时，获取GNSS网序贯平差算法中，前后期平差模型的 共同的参数和独立的参数，包括：前期子网独立站点坐标参数、子网间公共站点 坐标参数以及后期子网新增站点坐标参数；构建前后期平差模型的第一组误差方 程和第二组误差方程；In the specific implementation of this embodiment, in the sequential adjustment algorithm of the GNSS network, the common parameters and independent parameters of the previous and later adjustment models are obtained, including: the coordinate parameters of the independent sites of the sub-networks in the early stage, the coordinate parameters of the public sites between the sub-networks, and Added site coordinate parameters for later sub-networks; builds the first set of error equations and the second set of error equations for the adjustment model before and after;

对所述第一组误差方程进行单独平差，获取得到参数改正值矩阵和参数协方 差矩阵；Perform a separate adjustment on the first group of error equations to obtain a parameter correction value matrix and a parameter covariance matrix;

根据所述参数改正值矩阵计算得到未知参数的第一次平差值；Calculate the first adjustment value of the unknown parameter according to the parameter correction value matrix;

对参数协方差阵取对角值计算得到公共参数的中误差；Take the diagonal value of the parameter covariance matrix to calculate the median error of the common parameter;

将所述第一次平差值作为第二次平差时的近似值，代入所述第二组误差方程 中，计算新的常数项，定义新的观测信息改正数为V′2，得到新的误差方程；Taking the first adjustment value as the approximate value of the second adjustment, and substituting it into the second set of error equations, calculating a new constant term, defining the new observation information correction number as V′ 2 , and obtaining a new error equation;

根据模糊理论，以所述公共参数第一次平差值和所述中误差确定参数的模糊 中心和模糊幅度，根据所述常数项、所述模糊中心和所述模糊幅度构建平差函数 约束模型；According to fuzzy theory, the first adjustment value of the common parameter and the intermediate error determine the fuzzy center and the fuzzy amplitude of the parameter, and the adjustment function constraint model is constructed according to the constant term, the fuzzy center and the fuzzy amplitude ;

对所述构建平差约束函数模型求解得到参数的第二改正值；The second correction value of the parameter is obtained by solving the built adjustment constraint function model;

根据所述第二改正值和所述第一次平差值计算得到未知参数的第二次平差 值；Calculate the second adjustment value of the unknown parameter according to the second correction value and the first adjustment value;

根据所述第二次平差值计算获得各GNSS网点的坐标。The coordinates of each GNSS network point are obtained by calculating according to the second adjustment value.

通过对序贯平差加以改进，将前期平差得到的参数信息以约束条件的形式纳 入后期平差模型中解算，利用前期得到的先验信息对参数加以约束，提高模型的 抗误差干扰性。By improving the sequential adjustment, the parameter information obtained from the previous adjustment is incorporated into the later adjustment model in the form of constraints for calculation, and the prior information obtained in the early stage is used to constrain the parameters to improve the anti-error interference of the model. .

在本发明提供的又一实施例中，所述第一组误差方程为V1＝A11Xa+A12Xb-f1；In another embodiment provided by the present invention, the first set of error equations is V 1 =A 11 X a +A 12 X b -f 1 ;

所述第二组误差方程为V2＝B22Xb+B23Y-f2；The second set of error equations is V 2 =B 22 X b +B 23 Yf 2 ;

其中，V1为第一次观测值改正数，A11和A12为第一次平差模型系数矩阵，V2为第二次观测值改正数，B22和B23为第二次平差模型系数矩阵，Xa为所述前期子 网独立站点坐标参数，Xb为所述子网间公共站点坐标参数，Y为所述后期子网新 增站点坐标参数，f1为第一次平差误差方程的常数项，

Figure BDA0003334540950000081
f2为第二次平差误差方程的常数项，
Figure BDA0003334540950000082
L1和L2分别为第一 次平差观测值和第二次平差观测值，
Figure BDA0003334540950000083
和Y0为未知参数第一次平差解算时 所取的近似值。Among them, V 1 is the correction number of the first observation, A 11 and A 12 are the coefficient matrix of the first adjustment model, V 2 is the correction number of the second observation, and B 22 and B 23 are the second adjustment. Model coefficient matrix, X a is the coordinate parameter of the independent site of the sub-network in the early stage, X b is the coordinate parameter of the public site between the sub-networks, Y is the coordinate parameter of the newly added site of the sub-network in the later stage, and f 1 is the first flattening parameter. the constant term of the difference error equation,
Figure BDA0003334540950000081
f 2 is the constant term of the second adjustment error equation,
Figure BDA0003334540950000082
L 1 and L 2 are the first adjustment observations and the second adjustment observations, respectively,
Figure BDA0003334540950000083
and Y 0 are the approximations taken during the first adjustment of the unknown parameters.

在本实施例具体实施时，通过GNSS网序贯平差算法中，构建第一组误差方 程V1和第二组误差方程V2；In the specific implementation of this embodiment, the first group of error equations V 1 and the second group of error equations V 2 are constructed through the GNSS network sequential adjustment algorithm;

其中，V1＝A11Xa+A12Xb-f1，V2＝B22Xb+B23Y-f2，V1为第一次观 测值改正数，A11、A12为第一次平差模型系数阵，f1为其常数项

Figure BDA0003334540950000084
V2为第二次观测值改正数，B22、B23为第二次平差模 型系数阵，f2为其常数项
Figure BDA0003334540950000085
Xa为前期子网独立站点坐标参 数，Xb为子网间公共站点坐标参数，Y为后期子网新增站点坐标参数；L1、L2为 第一次、第二次平差观测值，A11、A12、B22、B23分别为其系数阵，
Figure BDA0003334540950000086
Y0为 参数第一次参与平差解算时所取的近似值。Wherein, V 1 =A 11 X a +A 12 X b -f 1 , V 2 =B 22 X b +B 23 Yf 2 , V 1 is the correction number of the first observation, A 11 and A 12 are the first sub-adjustment model coefficient matrix, f 1 is its constant term
Figure BDA0003334540950000084
V 2 is the correction number of the second observation value, B 22 and B 23 are the coefficient matrix of the second adjustment model, and f 2 is the constant term
Figure BDA0003334540950000085
X a is the coordinate parameter of the independent site of the sub-network in the early stage, X b is the coordinate parameter of the common site between the sub-networks, and Y is the coordinate parameter of the newly-added site of the sub-network in the later stage; L 1 , L 2 are the first and second adjustment observations , A 11 , A 12 , B 22 , and B 23 are their coefficient matrices, respectively,
Figure BDA0003334540950000086
Y 0 is the approximate value taken by the parameter when it participates in the adjustment solution for the first time.

在本发明提供的又一实施例中，所述参数改正值矩阵为

Figure BDA0003334540950000087
In another embodiment provided by the present invention, the parameter correction value matrix is
Figure BDA0003334540950000087

所述参数协方差矩阵为

Figure BDA0003334540950000088
The parameter covariance matrix is
Figure BDA0003334540950000088

其中，

Figure BDA0003334540950000089
和
Figure BDA00033345409500000810
为未知参数的改正数，P1为观测值权矩阵，
Figure BDA00033345409500000811
为单位 权方差,r为第一次平差时多余观测数。。in,
Figure BDA0003334540950000089
and
Figure BDA00033345409500000810
is the correction number of the unknown parameter, P 1 is the observation weight matrix,
Figure BDA00033345409500000811
is the unit weight variance, and r is the number of excess observations in the first adjustment. .

在本实施例具体实施时，对第一组误差方程V1单独平差，获得参数改正值 矩阵

Figure BDA0003334540950000091
In the specific implementation of this embodiment, the first group of error equations V 1 are adjusted individually to obtain a parameter correction value matrix
Figure BDA0003334540950000091

并获得参数协方差矩阵为

Figure BDA0003334540950000092
and obtain the parameter covariance matrix as
Figure BDA0003334540950000092

式中，

Figure BDA0003334540950000093
为参数改正数，P1为观测值权阵，
Figure BDA0003334540950000094
为单位权方差,r为 第一次平差时多余观测数。In the formula,
Figure BDA0003334540950000093
is the parameter correction number, P 1 is the observation weight matrix,
Figure BDA0003334540950000094
is the unit weight variance, and r is the number of excess observations in the first adjustment.

在本发明提供的又一实施例中，所述第一次平差值为

Figure BDA0003334540950000095
In another embodiment provided by the present invention, the first adjustment value is
Figure BDA0003334540950000095

在本实施例具体实施时，根据所述参数改正值矩阵计算得到未知参数的第一 次平差值

Figure BDA0003334540950000096
During the specific implementation of this embodiment, the first adjustment value of the unknown parameter is calculated according to the parameter correction value matrix.
Figure BDA0003334540950000096

其中，

Figure BDA0003334540950000097
为第二次平差未知参数平差值，
Figure BDA0003334540950000098
为第一次平差未知参数的平差值，x″b，y′为第二次平差时参数的改正数，
Figure BDA0003334540950000099
为公共参数的权阵， 由第一次平差结果得到。B22、B23为第二次平差模型系数阵，P2为第二次观测值权 阵，l2为其常数项。in,
Figure BDA0003334540950000097
is the adjustment value of the unknown parameter for the second adjustment,
Figure BDA0003334540950000098
is the adjustment value of the unknown parameter in the first adjustment, x″ b , y′ is the correction number of the parameter in the second adjustment,
Figure BDA0003334540950000099
is the weight matrix of common parameters, obtained from the first adjustment result. B 22 and B 23 are the coefficient matrix of the second adjustment model, P 2 is the weight matrix of the second observation, and l 2 is the constant term.

在本发明提供的又一实施例中，所述协方差为

Figure BDA00033345409500000910
In yet another embodiment provided by the present invention, the covariance is
Figure BDA00033345409500000910

由协方差可知，传统序贯平差根据前期平差后的参数估值

Figure BDA00033345409500000911
及其协方差阵
Figure BDA00033345409500000912
结合当期观测数据进行整体平差，虽不需要前期观测值即可实现与整体平差一致 的解算效果。但若前期先验参数或当期观测信息含有粗差，则势必造成后及参数 及其协方差阵的扭曲。为削弱先验参数异常、观测粗差对参数估值的影响，本发 明结合GNSS网对序贯平差加以改进，将前期平差得到的参数信息以约束条件的 形式纳入后期平差模型中解算，利用前期得到的先验信息对参数加以约束，提高 模型的抗误差干扰性。It can be seen from the covariance that the traditional sequential adjustment is estimated based on the parameters after the previous adjustment.
Figure BDA00033345409500000911
and its covariance matrix
Figure BDA00033345409500000912
Combined with the current observation data to carry out the overall adjustment, although the previous observation value is not required, the solution effect consistent with the overall adjustment can be achieved. However, if the prior prior parameters or current observation information contains gross errors, it will inevitably lead to distortion of the posterior parameters and their covariance matrix. In order to weaken the influence of abnormal a priori parameters and gross observation errors on parameter estimation, the present invention improves the sequential adjustment in combination with the GNSS network, and incorporates the parameter information obtained from the previous adjustment into the later adjustment model in the form of constraints to solve the problem. It uses the prior information obtained in the early stage to constrain the parameters to improve the anti-error interference of the model.

在本发明提供的又一实施例中，所述将所述第一次平差值作为第二次平差时 的近似值，代入所述第二组误差方程中，计算新的常数项，得到新的误差方程， 具体包括：In another embodiment provided by the present invention, the first adjustment value is used as an approximate value during the second adjustment, and is substituted into the second set of error equations to calculate a new constant term to obtain a new The error equation of , specifically includes:

将所述第一次平差值

Figure BDA0003334540950000101
作为第二次平差时的近似值，代入所述第二组平差 误差方程中计算得到新的常数项l2，定义新的观测信息改正数为V′2，得到新的误 差方程；the first adjustment value
Figure BDA0003334540950000101
As an approximate value during the second adjustment, a new constant term l 2 is calculated by substituting it into the second set of adjustment error equations, and a new correction number of observation information is defined as V′ 2 to obtain a new error equation;

其中，

Figure BDA0003334540950000102
in,
Figure BDA0003334540950000102

在本实施例具体实施时，In the specific implementation of this embodiment,

根据前期平差结果，我们可以认为测站坐标(即未知参数)处于一定的空间 范围内，可以形象的表示为以前期的参数估值为中心，以Δ为半径的范围内，本 结合模糊理论，将前期参数平差得到的信息构造为一个模糊数：X∈[X0 Δ]， X～μ(X)；According to the previous adjustment results, we can think that the coordinates of the station (that is, the unknown parameters) are within a certain spatial range, which can be visually expressed as the center of the previous parameter estimates and the range with Δ as the radius. This combination of fuzzy theory , construct the information obtained from the previous parameter adjustment as a fuzzy number: X∈[X 0 Δ], X～μ(X);

X0为模糊中心，可取为参数的前期平差估值，Δ为模糊范围，参数可取 Δ＝n*σ，σ为参数中误差，可由前期平差结果中得到。在测量误差理论中，偶 然误差服从正态分布，且P{-3σ<Δ<3σ}＝99.7％，即误差超过3倍中误差的概率 为0.3％，是小概率事件，因此可取n＝3。μ(X)为隶属函数，表示元素隶属于模糊 集合的程度，常见的隶属函数有正态模糊数、正弦模糊数等。X 0 is the fuzzy center, which can be taken as the estimation of the previous adjustment of the parameter, Δ is the fuzzy range, the parameter can be taken as Δ=n*σ, and σ is the error in the parameter, which can be obtained from the previous adjustment result. In the measurement error theory, the accidental error obeys a normal distribution, and P{-3σ<Δ<3σ}=99.7%, that is, the probability of the error exceeding 3 times is 0.3%, which is a small probability event, so n=3 . μ(X) is a membership function, which indicates the degree to which an element belongs to a fuzzy set. Common membership functions include normal fuzzy numbers, sinusoidal fuzzy numbers, and so on.

在本发明提供的又一实施例中，所述根据模糊理论，以所述第一次平差值和 所述中误差确定参数的模糊中心和模糊幅度，根据所述常数项、所述模糊中心和 所述模糊幅度构建平差函数约束模型，具体包括：In another embodiment provided by the present invention, according to the fuzzy theory, the fuzzy center and the fuzzy magnitude of the parameters are determined by the first adjustment value and the intermediate error, and the fuzzy center and the fuzzy center are determined according to the constant term, the fuzzy center and the fuzzy magnitude to construct an adjustment function constraint model, which specifically includes:

以公共参数的第一次平差值

Figure BDA0003334540950000103
作为参数的模糊中心，则参数改正值的模糊中 心
Figure BDA0003334540950000111
为第二次平差时参数所取近似值。以所述中误差的3倍值为 模糊幅度Δ前；Take the first adjustment value of the common parameter
Figure BDA0003334540950000103
As the fuzzy center of the parameter, the fuzzy center of the parameter correction value
Figure BDA0003334540950000111
Approximate values for the parameters in the second adjustment. Take 3 times of the middle error as the blur amplitude Δ before ;

根据隶属函数、所述模糊中心和所述模糊幅度构建所述平差函数模型：

Figure BDA0003334540950000112
The adjustment function model is constructed according to the membership function, the fuzzy center and the fuzzy magnitude:
Figure BDA0003334540950000112

其中，x″b和y′为第二次平差时参数的改正数，μA(x″b)为x″b的隶属函数，

Figure BDA0003334540950000113
，V′2为第二次平差观测信息改正数，B22、B23为第二次 平差模型系数阵，l2为其常数项。Among them, x″ b and y′ are the correction numbers of the parameters during the second adjustment, μ A (x″ b ) is the membership function of x″ b ,
Figure BDA0003334540950000113
, V′ 2 is the correction number of the second adjustment observation information, B 22 and B 23 are the second adjustment model coefficient matrix, and l 2 is the constant term.

平差函数模型可理解为部分参数带有约束条件的平差模型，μ(x″b)为隶属函 数，表示元素隶属于模糊数的程度，以正态模糊数为例，参数的可能性分布可表 示为

Figure BDA0003334540950000114
The adjustment function model can be understood as an adjustment model with some parameters with constraints, μ(x″ b ) is the membership function, which indicates the degree of membership of the element to the fuzzy number. Taking the normal fuzzy number as an example, the probability distribution of the parameter can be expressed as
Figure BDA0003334540950000114

在本发明提供的又一实施例中，所述对所述构建平差约束函数模型求解得到 参数的第二改正值，具体包括：In another embodiment provided by the present invention, the second correction value of the parameter obtained by solving the described construction adjustment constraint function model specifically includes:

对观测残差平方和取最小值的同时，x″b的隶属函数μA(x″b)取最大值，得到 准则函数

Figure BDA0003334540950000115
While taking the minimum value of the sum of squares of the observed residuals, the membership function μ A (x″ b ) of x″ b takes the maximum value to obtain the criterion function
Figure BDA0003334540950000115

根据所述模糊幅度建立算子

Figure BDA0003334540950000116
Build an operator based on the blurring magnitude
Figure BDA0003334540950000116

根据所述算子将所述准则函数转化为准则函数矩阵

Figure BDA0003334540950000117
Transform the criterion function into a criterion function matrix according to the operator
Figure BDA0003334540950000117

对所述准则函数矩阵求偏导并其等于0，计算得到参数的第二改正值

Figure BDA0003334540950000118
Calculate the partial derivative of the criterion function matrix and make it equal to 0, and calculate the second corrected value of the parameter
Figure BDA0003334540950000118

其中，τ为任一数值，0<τ<1，W＝diag[w1 w2 … wt]，Pi为观测值权阵，

Figure BDA0003334540950000121
为观测值残差，n＝1，2，3…，t＝1，2，3…，j＝1，2，...，t，Vxb＝x″b-xb前，表 示参数改正值与其先验模糊中心的偏差。Among them, τ is any value, 0<τ<1, W=diag[w 1 w 2 … w t ], P i is the observation weight matrix,
Figure BDA0003334540950000121
is the residual error of the observation value, n=1, 2, 3..., t=1, 2, 3..., j=1, 2,..., t, V xb = x″ b -x b before , it means the parameter correction The deviation of the value from its prior fuzzy center.

在本实施例具体实施时，构建算子

Figure BDA0003334540950000122
为避免出现Δ前＝0导致wj无限大 的情况，可将算子设定为
Figure BDA0003334540950000123
0<τ<1，为一适当小的数。During the specific implementation of this embodiment, the operator is constructed
Figure BDA0003334540950000122
In order to avoid the situation where w j is infinitely large due to Δbefore = 0, the operator can be set as
Figure BDA0003334540950000123
0<τ<1, is an appropriately small number.

取W＝diag[w1 w2 … wt]；Take W=diag[w 1 w 2 ... w t ];

则准则函数可写为准则函数矩阵

Figure BDA0003334540950000124
Then the criterion function can be written as the criterion function matrix
Figure BDA0003334540950000124

对准则函数矩阵求偏导并令其等于零可求出参数解；The parametric solution can be obtained by taking the partial derivative of the criterion function matrix and making it equal to zero;

参数解

Figure BDA0003334540950000125
parametric solution
Figure BDA0003334540950000125

其中，x″b，y′为第二次平差时参数的改正数，P2为第二次平差观测值权阵， B22、B23为系数阵，l2为常数项，xb前为参数改正数模糊中心。Among them, x″ b , y′ are the correction numbers of the parameters during the second adjustment, P 2 is the weight matrix of the second adjustment observations, B 22 and B 23 are the coefficient matrices, l 2 is the constant term, x b The former is the parameter correction number fuzzy center.

在本发明提供的又一实施例中，所述根据所述第二改正值和所述第一次平差 值计算得到未知参数的第二次平差值，具体为：In another embodiment provided by the present invention, the second adjustment value of the unknown parameter is calculated according to the second correction value and the first adjustment value, specifically:

将所述第二改正值

Figure BDA0003334540950000126
所述第一 次平差值
Figure BDA0003334540950000127
代入平差值计算公式计算第二次平差值；the second correction value
Figure BDA0003334540950000126
the first adjustment value
Figure BDA0003334540950000127
Substitute into the adjustment value calculation formula to calculate the second adjustment value;

所述平差值计算公式为

Figure BDA0003334540950000128
The calculation formula of the adjustment value is:
Figure BDA0003334540950000128

在本实施例具体实施时，将所述第一次平差值

Figure BDA0003334540950000129
和参数解
Figure BDA00033345409500001210
代入平差值计算公式
Figure BDA0003334540950000131
计算第二次平差值
Figure BDA0003334540950000132
并根据计算的第二次平差值用于GNSS 基准站网解算，获得各GNSS网点的坐标。During the specific implementation of this embodiment, the first adjustment value is
Figure BDA0003334540950000129
and parametric solution
Figure BDA00033345409500001210
Substitute into the calculation formula of the adjustment value
Figure BDA0003334540950000131
Calculate the second adjustment value
Figure BDA0003334540950000132
According to the calculated second adjustment value, it is used for the calculation of the GNSS reference station network, and the coordinates of each GNSS network point are obtained.

在本发明提供的又一实施例中，将所述改进的GNSS网序贯平差计算方法应 用于GNSS网，参见图2所示，是本发明实施例提供的GNSS网的网型图；In another embodiment provided by the present invention, the described improved GNSS network sequential adjustment calculation method is applied to the GNSS network, referring to shown in Figure 2, it is the network diagram of the GNSS network provided by the embodiment of the present invention;

采用两台GNSS接收机进行同步观测，LC01、LC03为已知点的接收机，LC02、 LC04视为未知点进行平差计算；Two GNSS receivers are used for synchronous observation, LC01 and LC03 are receivers with known points, and LC02 and LC04 are regarded as unknown points for adjustment calculation;

LC01、LC03、LC02、LC04四个测站的坐标理论值如表1：The theoretical coordinate values of the four stations LC01, LC03, LC02 and LC04 are shown in Table 1:

表1四个测站坐标理论值Table 1 Theoretical values of the coordinates of the four stations

Figure BDA0003334540950000133
Figure BDA0003334540950000133

第一期选择1、2、3三个基线向量，第二期选择4、5两个基线向量。通过 MATLAB仿真系统对各基线坐标差理论值附加偶然误差并对基线4、5附加粗差， 形成观测基线信息如表2所示。解算过程中基线方差阵取单位矩阵。Three baseline vectors 1, 2, and 3 were selected in the first phase, and two baseline vectors 4 and 5 were selected in the second phase. Through the MATLAB simulation system, random errors are added to the theoretical value of each baseline coordinate difference, and gross errors are added to the baselines 4 and 5 to form the observation baseline information as shown in Table 2. During the solution process, the baseline variance matrix takes the identity matrix.

表2观测基线信息Table 2 Observation baseline information

Figure BDA0003334540950000134
Figure BDA0003334540950000134

分别利用序贯最小二乘及约束序贯算法进行解算，坐标解算结果及坐标残差 分别如表3、表4所示；Sequential least squares and constrained sequential algorithms are used to solve respectively, and the coordinate solution results and coordinate residuals are shown in Table 3 and Table 4 respectively;

表3未知点坐标解算结果/mTable 3 Unknown point coordinate solution results/m

Figure BDA0003334540950000135
Figure BDA0003334540950000135

表4坐标残差/mTable 4 Coordinate residuals/m

Figure BDA0003334540950000141
Figure BDA0003334540950000141

并将序贯最小二乘及约束序贯算法的坐标误差和坐标残差分别进行对比，如 表5和图3所示；The coordinate errors and coordinate residuals of the sequential least squares and constrained sequential algorithms are compared respectively, as shown in Table 5 and Figure 3;

表5未知点坐标误差Table 5 Coordinate error of unknown point

Figure RE-GDA0003589279770000142
Figure RE-GDA0003589279770000142

通过分析以上结果，可以得到：约束序贯解法的坐标残差分别为0.0215m及0.0229m，均小于序贯最小二乘的解算结果；并且约束序贯解法得到的的坐标估 值与理论值的残差均小于序贯最小二乘。说明约束序贯解法充分利用了前期平差 得到的参数先验信息，较于传统最小二乘，具有较好的抗粗差干扰性。By analyzing the above results, it can be obtained that the coordinate residuals of the constrained sequential solution method are 0.0215m and 0.0229m respectively, which are both smaller than the solution results of the sequential least squares; The residuals are smaller than sequential least squares. It shows that the constrained sequential solution method makes full use of the prior information of the parameters obtained by the previous adjustment, and has better resistance to gross error interference than the traditional least squares method.

本发明提出的约束序贯平差模型在保证解算效率的同时提高GNSS网的抗误 差干扰性，提升GNSS网点坐标的解算精度。可应用于以下领域：大规模GNSS 基准站网解算；GNSS技术进行控制网分期布设数据处理。The constrained sequential adjustment model proposed by the present invention improves the anti-error interference of the GNSS network while ensuring the calculation efficiency, and improves the calculation accuracy of the coordinates of the GNSS network point. It can be applied to the following fields: large-scale GNSS reference station network solution; GNSS technology for control network staged data processing.

本发明提供一种改进的GNSS网序贯平差计算方法，通过建立前后期平差模 型的第一组误差方程和第二组误差方程；对所述第一组误差方程单独进行平差， 得到参数改正值矩阵和参数协方差矩阵；根据所述参数改正值矩阵计算得到未知 参数的第一次平差值；对参数协方差阵取对角值计算得到公共参数的中误差；根 据模糊理论，以所述公共参数第一次平差值和所述中误差确定参数的模糊中心和 模糊幅度，根据所述常数项、所述模糊中心和所述模糊幅度构建平差函数约束模 型；对所述构建平差约束函数模型求解得到参数的第二改正值；根据所述第二改 正值和所述第一次平差值计算得到未知参数的第二次平差值；根据所述第二次平 差值计算获得各GNSS网点的坐标。当后期观测信息含有粗差时，可以有效削弱 粗差带来的参数估值扭曲，减少误差积累，提高解算精度。The present invention provides an improved GNSS network sequential adjustment calculation method. By establishing the first group of error equations and the second group of error equations of the adjustment models before and after the adjustment, and adjusting the first group of error equations separately, the result is obtained: The parameter correction value matrix and the parameter covariance matrix; the first adjustment value of the unknown parameter is calculated according to the parameter correction value matrix; the diagonal value of the parameter covariance matrix is calculated to obtain the intermediate error of the public parameter; according to the fuzzy theory, Determine the fuzzy center and fuzzy amplitude of the parameter with the first adjustment value of the common parameter and the intermediate error, and construct an adjustment function constraint model according to the constant term, the fuzzy center and the fuzzy amplitude; Building an adjustment constraint function model and solving to obtain the second correction value of the parameter; calculating the second adjustment value of the unknown parameter according to the second correction value and the first adjustment value; according to the second adjustment value The difference calculation obtains the coordinates of each GNSS network point. When the later observation information contains gross errors, the distortion of parameter estimation caused by gross errors can be effectively weakened, the accumulation of errors can be reduced, and the calculation accuracy can be improved.

以上所述是本发明的优选实施方式，应当指出，对于本技术领域的普通技术 人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改 进和润饰也视为本发明的保护范围。The above are the preferred embodiments of the present invention. It should be pointed out that for those skilled in the art, without departing from the principles of the present invention, several improvements and modifications can also be made, and these improvements and modifications may also be regarded as It is the protection scope of the present invention.

Claims (9)

1. An improved GNSS network sequential adjustment calculation method is characterized by comprising the following steps:
establishing a first group of error equations and a second group of error equations of a front-stage adjustment model and a rear-stage adjustment model according to the geometrical relationship of a baseline vector among measuring stations by defining coordinate parameters of an early-stage subnet independent measuring station, coordinate parameters of an inter-subnet common measuring station and coordinate parameters of a later-stage subnet independent measuring station in the GNSS network;
performing adjustment on the first group of error equations separately to obtain a parameter correction value matrix and a parameter covariance matrix;
calculating according to the parameter correction value matrix to obtain a first adjustment value of the unknown parameter;
calculating a diagonal value of the parameter covariance matrix to obtain a medium error of an unknown parameter;
substituting the first adjustment value of the public parameter as an approximate value of a second adjustment value into the second adjustment error equation, calculating a new constant term, redefining an observed value correction number, and obtaining a new error equation;
determining a fuzzy center and a fuzzy amplitude of a correction value of a common parameter according to the fuzzy theory and the first adjustment value and the median error, and constructing an adjustment function constraint model according to the constant term, the fuzzy center and the fuzzy amplitude;
solving the constructed adjustment constraint function model to obtain a second correction value of the parameter;
calculating a second adjustment value of the unknown parameter according to the second correction value and the first adjustment value;
and calculating to obtain the coordinates of each GNSS network point according to the second adjustment value.
2. The method for improved GNSS network sequential adjustment computation of claim 1, wherein the first set of adjustment error equations is V1＝A11Xa+A12Xb-f1；
The second set of adjustment error equations is V2＝B22Xb+B23Y-f2；
Wherein, V1For first observation correction, A11And A12Is a first adjustment model coefficient matrix, V2Number of second observation correction, B22And B23Is a second adjustment model coefficient matrix, XaIs the coordinate parameter, X, of the preceding subnet independent stationbFor the coordinate parameter of the common station between the subnets, Y is the coordinate parameter of the newly added station of the later subnet, f1Is a constant term of the first adjustment error equation,
Figure FDA0003334540940000021
f2is a constant term of the second adjustment error equation,
Figure FDA0003334540940000022
L1and L2A first set of observations and a second set of observations,
Figure FDA0003334540940000023
and Y0For the first time for unknown parametersAnd (4) an approximate value taken in the adjustment calculation.
3. The improved GNSS network sequential adjustment calculation method of claim 2, wherein the parameter correction value matrix is
Figure FDA0003334540940000024
The parameter covariance matrix is
Figure FDA0003334540940000025
Wherein,
Figure FDA0003334540940000026
and
Figure FDA0003334540940000027
as a correction of an unknown parameter, P1In order to obtain a matrix of weights of the observations,
Figure FDA0003334540940000028
is the variance of the unit weight, and r is the number of redundant observations at the first adjustment.
4. The method as recited in claim 3, wherein the first adjustment value is the sequential adjustment of the GNSS network
Figure FDA0003334540940000029
5. The method for improved GNSS network sequential adjustment computation of claim 4, wherein the covariance is
Figure FDA00033345409400000210
6. The improved GNSS network sequential adjustment calculation method according to claim 5, wherein the substituting the first adjustment value of the common parameter as an approximate value of a second set of adjustments into the second set of adjustment error equations to calculate a new constant term, redefining an observation value correction, and obtaining a new error equation specifically includes:
the adjustment value of the common parameter in the first adjustment
Figure FDA0003334540940000031
As an approximate value in the second adjustment, substituting the approximate value into the second adjustment error equation, and calculating a new constant term l2Defining a new observation information correction number as V2', obtaining a new error equation;
wherein,
Figure FDA0003334540940000032
7. the improved GNSS network sequential adjustment calculation method according to claim 6, wherein the determining a fuzzy center and a fuzzy magnitude of parameter correction values according to the fuzzy theory with the first adjustment value and the median error, and constructing an adjustment function constraint model according to the constant term, the fuzzy center and the fuzzy magnitude specifically comprises:
by first adjustment of a common parameter
Figure FDA0003334540940000033
As the fuzzy center of the parameter, the fuzzy center of the parameter correction value
Figure FDA0003334540940000034
Figure FDA0003334540940000035
And the parameters are taken as approximate values during the second adjustment. Taking the value of 3 times of the medium error as the fuzzy amplitude deltaFront side；
Constructing the constraint according to membership functions, the fuzzy center and the fuzzy amplitudeAdjustment function model:
Figure FDA0003334540940000036
wherein, x ″)bAnd y' is the correction of the parameter at the second adjustment, μA(x″b) Is x ″)bIs a function of the membership of (a) to (b),
Figure FDA0003334540940000037
8. the improved GNSS network sequential adjustment calculation method according to claim 6, wherein the solving the constructed constrained adjustment function model to obtain a second correction value of the parameter specifically includes:
taking the minimum value of the sum of squares of the observed residuals, x ″)bMembership function μ ofA(x″b) Taking the maximum value to obtain a criterion function
Figure FDA0003334540940000038
Establishing an operator from the fuzzy magnitude
Figure FDA0003334540940000041
Converting the criterion function into a criterion function matrix according to the operator
Figure FDA0003334540940000042
Calculating the deviation of the criterion function matrix and making it equal to 0, calculating to obtain the second correction value of the parameter
Figure FDA0003334540940000043
Wherein τ is any number, 0<τ<1，W＝diag[w1 w2 … wt]，PiIn order to obtain a weighted array of observations,
Figure FDA0003334540940000044
for the observed residuals, n is 1, 2, 3 …, t is 1, 2, 3 …, j is 1, 2xb＝x″b-xb front ofAnd represents the deviation of the parameter correction value from its a priori blur center.
9. The improved GNSS network sequential adjustment calculation method according to claim 8, wherein the second adjustment value of the unknown parameter is calculated according to the second correction value and the first adjustment value, and specifically includes:
correcting the second value
Figure FDA0003334540940000045
And the first adjustment value
Figure FDA0003334540940000046
Substituting the average value into an average value calculation formula to calculate a secondary average value;
the average value is calculated by the formula
Figure FDA0003334540940000047
CN202111290536.8A 2021-11-02 2021-11-02 Improved GNSS network sequential adjustment calculation method Active CN114662268B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111290536.8A CN114662268B (en) 2021-11-02 2021-11-02 Improved GNSS network sequential adjustment calculation method

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111290536.8A CN114662268B (en) 2021-11-02 2021-11-02 Improved GNSS network sequential adjustment calculation method

Publications (2)

Publication Number Publication Date
CN114662268A true CN114662268A (en) 2022-06-24
CN114662268B CN114662268B (en) 2023-04-07

Family

ID=82026227

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111290536.8A Active CN114662268B (en) 2021-11-02 2021-11-02 Improved GNSS network sequential adjustment calculation method

Country Status (1)

Country Link
CN (1) CN114662268B (en)

Citations (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20130342393A1 (en) * 2012-06-20 2013-12-26 Topcon Positioning Systems, Inc. Selection of a Subset of Global Navigation Satellite System Measurements Based on Relation between Shifts in Target Parameters and Sum of Residuals
CN108802781A (en) * 2018-05-03 2018-11-13 广州市中海达测绘仪器有限公司 The acquisition methods of integer ambiguity
CN109061641A (en) * 2018-07-06 2018-12-21 中南大学 A kind of InSAR timing earth's surface deformation monitoring method based on sequential adjustment
CN109085629A (en) * 2018-08-29 2018-12-25 广州市中海达测绘仪器有限公司 GNSS baseline vector procession localization method, device and navigation equipment
CN109633723A (en) * 2018-12-26 2019-04-16 东南大学 A kind of single epoch GNSS calculation method of attached horizontal restraint
CN110673182A (en) * 2019-09-29 2020-01-10 清华大学 A kind of GNSS high-precision and fast positioning method and device
CN113295149A (en) * 2021-05-17 2021-08-24 中铁第四勘察设计院集团有限公司 CP III coordinate calculation method and device based on joint observation quantity
CN113325453A (en) * 2021-06-22 2021-08-31 中国科学院精密测量科学与技术创新研究院 GNSS non-differential ambiguity determination method based on parameter constraint and rapid positioning method
CN113343163A (en) * 2021-04-19 2021-09-03 华南农业大学 Large-scale corner mesh adjustment and precision evaluation method, system and storage medium
CN113358017A (en) * 2021-06-02 2021-09-07 同济大学 Multi-station cooperative processing GNSS high-precision deformation monitoring method

Patent Citations (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20130342393A1 (en) * 2012-06-20 2013-12-26 Topcon Positioning Systems, Inc. Selection of a Subset of Global Navigation Satellite System Measurements Based on Relation between Shifts in Target Parameters and Sum of Residuals
CN108802781A (en) * 2018-05-03 2018-11-13 广州市中海达测绘仪器有限公司 The acquisition methods of integer ambiguity
CN109061641A (en) * 2018-07-06 2018-12-21 中南大学 A kind of InSAR timing earth's surface deformation monitoring method based on sequential adjustment
CN109085629A (en) * 2018-08-29 2018-12-25 广州市中海达测绘仪器有限公司 GNSS baseline vector procession localization method, device and navigation equipment
CN109633723A (en) * 2018-12-26 2019-04-16 东南大学 A kind of single epoch GNSS calculation method of attached horizontal restraint
CN110673182A (en) * 2019-09-29 2020-01-10 清华大学 A kind of GNSS high-precision and fast positioning method and device
CN113343163A (en) * 2021-04-19 2021-09-03 华南农业大学 Large-scale corner mesh adjustment and precision evaluation method, system and storage medium
CN113295149A (en) * 2021-05-17 2021-08-24 中铁第四勘察设计院集团有限公司 CP III coordinate calculation method and device based on joint observation quantity
CN113358017A (en) * 2021-06-02 2021-09-07 同济大学 Multi-station cooperative processing GNSS high-precision deformation monitoring method
CN113325453A (en) * 2021-06-22 2021-08-31 中国科学院精密测量科学与技术创新研究院 GNSS non-differential ambiguity determination method based on parameter constraint and rapid positioning method

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
刘洋 等: "组合GNSS系统定位数据质量与精度比较分析", 《北京测绘》 *
张明 等: "基于序贯平差的长距离基准站间模糊度快速固定", 《武汉大学学报(信息科学版)》 *
褚成凤 等: "基于方差-协方差阵的宝林隧洞控制点稳定性分析", 《人民长江》 *

Also Published As

Publication number Publication date
CN114662268B (en) 2023-04-07

Similar Documents

Publication Publication Date Title
WO2021203871A1 (en) Cooperative positioning method and apparatus, device, and storage medium
CN105629263B (en) A kind of troposphere atmosphere delay estimation error correcting method and correction system
EP4459327A1 (en) Terminal positioning method and apparatus, and device and medium
CN108985373B (en) Multi-sensor data weighting fusion method
CN111556454B (en) A Weighted DV_Hop Node Location Method Based on Minimum Mean Square Error Criterion
CN110618438B (en) Atmospheric error calculation method and device, computer equipment and storage medium
CN108960334B (en) Multi-sensor data weighting fusion method
CN109782269B (en) Distributed multi-platform cooperative active target tracking method
WO2019024895A1 (en) Real-time estimation and quality control method for gnss satellite clock difference
CN109754013B (en) Electric power system hybrid measurement fusion method based on unscented Kalman filtering
CN110426717B (en) A cooperative positioning method and system, positioning device, and storage medium
CN112597428A (en) Flutter detection correction method based on beam adjustment and image resampling of RFM model
CN110493869B (en) RSSI-based K-nearest neighbor differential correction centroid positioning method
CN106353722A (en) RSSI (received signal strength indicator) distance measuring method based on cost-reference particle filter
CN112153564B (en) Efficient multi-hop positioning method based on combination of centralized and distributed computing
CN113295149A (en) CP III coordinate calculation method and device based on joint observation quantity
CN112865096A (en) Power distribution network state estimation method and system considering PMU (phasor measurement Unit) measurement phase angle deviation
CN108684074A (en) Distance measuring method based on RSSI and device
CN114662268A (en) An Improved GNSS Network Sequential Adjustment Calculation Method
CN111354040A (en) An adjustment method for optical satellite image block network based on Partial EIV model
CN112020005B (en) Method, device and system for solving Long Term Evolution (LTE) network mode three-interference
CN108848447B (en) A Differential DV_Distance Node Localization Method Using Unknown Node Correction
CN107909606A (en) A kind of SAR image registration communication center elimination of rough difference method
CN112835079A (en) GNSS self-adaptive weighting positioning method based on edge sampling consistency
CN110673088A (en) Time-of-arrival-based object localization in mixed line-of-sight and non-line-of-sight environments

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