WO2023019613A1

WO2023019613A1 - Large-scale gnss network parallel resolution method and system based on dynamic partioning

Info

Publication number: WO2023019613A1
Application number: PCT/CN2021/114428
Authority: WO
Inventors: 赵祥伟; 陈正宇; 李欣; 秦臻; 徐君民; 王宁; 陈杰; 杨晓慧
Original assignee: 中国能源建设集团江苏省电力设计院有限公司
Priority date: 2021-08-19
Filing date: 2021-08-25
Publication date: 2023-02-23
Also published as: ZA202307199B; CN113848577A

Abstract

A dynamic-partitioning-based large-scale GNSS network automatic parallel resolving method is applicable to a large-scale GNSS network resolving work field. The method comprises the steps of: 1) data collection and pre-processing; 2) automatic partitioning of an observation station subnet; 3) parallel automatic baseline resolution of partitions; and 4) comprehensive adjustment. The partitioning of the observation station and the parallel resolution of the partitions improve calculation efficiency.

Description

A large-scale GNSS network parallel computing method and system based on dynamic partition

technical field

The invention belongs to the technical field of satellite positioning, and in particular relates to a large-scale GNSS network parallel computing method based on dynamic partitioning, and also relates to a large-scale GNSS network parallel computing system based on dynamic partitioning.

Background technique

With the acceleration of the global coverage of my country's Beidou system and the encryption and upgrade of the ground-based reference station network, large-scale GNSS networks have ushered in new opportunities for development. However, while the GNSS original observation information is continuously enriched, it also proposes a series of problems for relevant professionals in the aspects of large-scale data preprocessing, organization management, reference station selection, data processing strategy, partitioning scheme, quality inspection and evaluation, etc. challenge. At the same time, in order to meet the needs of high-precision data processing, it is usually necessary to use special data processing software (such as GAMIT, GIPSY and Bernese) for calculation. At present, the accuracy of the single-day coordinate solution of the reference station reaches 3-5mm in the horizontal direction and 6-8mm in the elevation direction.

Limited by computer software and hardware, most GNSS data processing software can only process data from less than 100 stations at the same time (such as GAMIT). If too many stations are processed at the same time, a large amount of computer hardware resources and Time will seriously affect the efficiency of data calculation and lead to the lag of calculation results. The strategy usually adopted at present is to divide a large-scale GNSS reference station network into several subnetworks, each subnetwork is solved independently and then jointly processed to obtain the final result. However, when subnets are used to solve large-scale GNSS reference station network data, how to divide subnets and select public stations has a certain impact on the calculation results.

In recent years, with the advancement of network technology, new distributed computing technologies such as grid computing and cloud computing have emerged. a new space for development. Distributed computing technology is being gradually applied to the field of geodesy due to its efficient resource utilization and computing efficiency. Since the essence of high-performance computing technology is parallel computing, decomposing the traditional GNSS data processing technology and designing a reasonable parallel algorithm is the primary problem to realize the rapid processing of GNSS data.

Contents of the invention

The purpose of the present invention is to overcome the deficiencies in the prior art, and provide a large-scale GNSS network parallel computing method based on dynamic partitioning, which partitions observation stations and performs parallel computing for each partition, thereby improving computing efficiency.

In order to solve the problems of the technologies described above, the technical solutions provided by the invention are as follows:

In the first aspect, the present invention provides a large-scale GNSS network parallel solution method based on dynamic partition, including:

collect observational data;

Carry out subnetwork partitioning for the observation stations corresponding to the observation data to obtain a list of station partitions;

Based on the station partition list and observation data, parallel baseline calculations are performed for each partition to obtain a single-day baseline solution for each partition;

Adjustment processing is carried out on the single-day baseline solution of each division to obtain the final positioning results of each observation station.

Optionally, after the collection of observation data, it also includes

The observation data is preprocessed to obtain complete and qualified observation data products. The preprocessing includes data continuity check, integrity check, data editing and data quality check processing.

Optionally, the subnet partitioning method includes any one of geographic region partitioning, central diffusion partitioning, and central uniform partitioning.

Optionally, the baseline solution for each partition is obtained to obtain a single-day baseline solution for each partition, including:

Assuming that satellite i is observed at station j at time epoch _tj , the linearized dual-frequency carrier phase observation equation is written as:

Among them, the subscripts 1 and 2 represent the carrier phases of L1 and L2 respectively;

is the L1 carrier phase observation quantity, the unit is cycle;

is the L2 carrier phase observation quantity, the unit is cycle;

f ₁ is L1 carrier frequency;

f ₂ is the L2 carrier frequency;

τ _ij is the geometric propagation delay time of the carrier signal between the satellite and the receiver;

is the phase change caused by satellite and receiver clock differences;

is the tropospheric refraction delay on the propagation path of the carrier signal;

k _ij is the effect of ionospheric refraction;

v _ij is measurement error and residual error;

where n _ij is the integer ambiguity,

is the initial phase deviation of the receiver,

is the initial phase deviation of the satellite;

The geometric propagation delay is written as:

here

are the three-dimensional geocentric position vectors of the satellite and the station respectively, and c is the speed of light; the parameters to be solved are the geocentric position vectors of the station and the satellite;

If the phase change eliminated in the double-difference combination is omitted, the phase observation equation is simply written as:

In the formula, n _1ij is the integer number of L1, n _2ij is the integer number of L2;

If certain constraints are imposed on the effect of ionospheric refraction, an equation is also listed:

l _ij (t _j )＝k _ij (t _j )/f ₁ +v _kij (t _j ) (6)

In the formula, l _ij (t _j ) is the influence of ionospheric refraction at time t _j ; v _kij (t _j ) is v _1ij (t _j ) or v _2ij (t _j );

In each observation epoch, the observation equation of the above observations is written as the adjustment equation in matrix form:

In the formula,

is the true value of the observation,

is the coefficient matrix,

is the parameter to be estimated,

is the observation error correction number;

In the above formula are:

A _k =A _k1 =1/f ₁ ; A _k2 =1/f ₂ ; A _k1 =gA _k2 =A _k

The residual vector of double-differenced observations satisfies:

In the formula, D is a double-difference operator, which maps all phase one-way observations in the same epoch into independent double-difference observations, A is the linearized coefficient matrix, C is the cofactor matrix, and σ ² is the prior unit weight Variance, x _a is the unknown parameters to be solved, including the coordinate difference of observation stations, satellite orbit parameters, pole shift parameters, ambiguity parameters and tropospheric zenith delay parameters, x _k is the ionospheric delay parameters.

Optionally, after obtaining the single-day baseline solution for each partition, it also includes:

Perform a quality check on the baseline solution of each partition to obtain a qualified baseline result;

If there are unqualified conditions, check the quality of the original observation data. If it is caused by the poor quality of the station data, remove the station during the calculation.

Optionally, the specific process of the quality inspection is:

Calculate the evaluation index of the quality of the baseline results, the evaluation index is expressed by the standardized root mean square error NRMS, and the calculation formula is as follows:

In the formula, Y _i is the calculated value of each baseline vector epoch, Y is the true value of the baseline vector, N is the total number of epochs, n is the number of baselines,

is the variance of the calculated value for each epoch;

If the NRMS value is less than the set threshold, it is judged that the baseline solution is qualified, otherwise it is not qualified.

Optionally, the single-day baseline calculation of each subregion is adjusted to obtain the final positioning results of each observation station, including:

Let the position and velocity vectors of the observation station at time t ₀ be X ₀ ,

The position and velocity vectors of the station at time t _k are X _k ,

Then there are:

The third item on the right side of the above formula is the random drift part of the station between t ₀ and t _k , r is the random drift coefficient, δρ _k is the observation error at each moment, these influences are ignored, and the speed of the station is assumed does not change with time, formula (11) is written in matrix form as follows:

In the above formula, I is the unit matrix, and the formula (12) is considered as the state transition equation of the station position and velocity parameters, then the state transition matrix is:

Let the prior values of coordinates and velocity at this moment be X′ ₀ and

Then there are:

y _k is the parameter solution at time t _k , and A _k is the coefficient matrix at time t _k ;

The error equation is written as:

δx ₀ ,

is the correction amount of station position and velocity, V _k is obtained from the output result file;

The covariance matrix W _k of the random disturbance is a diagonal matrix. Since the drift of the station coordinates is subject to the Markov process, it is only necessary to give its power spectral density PSD value; the corresponding value is:

σ _i ² =PSD×(t _k+1 -t _k ) (17)

is the value of the i-th parameter in W _k , PSD is the empirical value;

Based on the constructed A _k , y _k , S _k and W _k , the Kalman filter estimation algorithm is used for adjustment processing, and the final positioning results of each observation station are obtained.

Optionally, after obtaining the final positioning results of each observation station, it also includes:

Perform quality checks on the final positioning results of each observation station to obtain the final qualified coordinate results of the GNSS station.

In the second aspect, the present invention also provides a large-scale GNSS network parallel computing system based on dynamic partitioning, including:

A data acquisition module for collecting observation data;

The partition processing module is used to carry out subnetwork partitions to the observation stations corresponding to the observation data, and obtain the station partition list;

The parallel calculation module is used to perform parallel baseline calculations for each zone based on the station zone list and observation data, and obtain a single-day baseline solution for each zone;

And, the comprehensive adjustment module is used for adjusting the single-day baseline calculation of each subregion to obtain the final positioning results of each observation station.

Compared with the prior art, the beneficial effects achieved by the present invention are: the present invention can automatically complete a series of tasks such as data download, subnet partition, baseline calculation, subnet fusion, and comprehensive adjustment without user interaction, thereby The calculation pressure is reduced, the calculation efficiency is improved, and the algorithm is not limited by the number of GNSS stations.

Description of drawings

In order to make the content of the present invention easier to understand clearly, the present invention will be described in further detail below according to specific embodiments in conjunction with the accompanying drawings, wherein:

Fig. 1 is the overall technical flow chart of the method of the present invention.

Detailed ways

The present invention will be further described below in conjunction with the accompanying drawings. The following examples are only used to illustrate the technical solution of the present invention more clearly, but not to limit the protection scope of the present invention.

In the description of the patent of the present invention, it should be noted that the terms "comprising", "comprising" or any other variant thereof are intended to cover a non-exclusive inclusion, except for those elements listed, and may also include elements not explicitly listed. other elements out.

Example 1

A kind of large-scale GNSS network parallel solution method based on dynamic partition of the present invention, flow process is shown in Fig. 1, comprises the following process:

The first step is to collect, download and organize the data such as GNSS observation data, broadcast ephemeris products, precise ephemeris data files (products released by the IGS Center), and calculation table files, and then preprocess the data (mainly refers to the data continuity check, completeness check, data editing and data quality check) to obtain complete and qualified data products.

Baseline calculation based on satellite observation data, satellite ephemeris and calculation table files to obtain high-precision baseline results. Due to the large number of measuring stations, it is necessary to perform partition calculation to improve the calculation efficiency, that is, to perform a separate baseline calculation process for the station data in each sub-region.

The second step is to extract the approximate coordinates of the observation station from the observation data file and convert them into latitude and longitude. Based on the name of the station and the latitude and longitude, the station is divided into subnetworks (partition methods include geographical area partition method, central diffusion partition method and central uniform partition method) Three partitioning schemes) to get the station name partition list (each partition has multiple stations).

This article mainly introduces three convenient and commonly used partitioning methods:

Geographical area division method: divide according to the location of the site's geographical location (commonly used mode, simple and fast);

Center diffusion partition method: according to the distance from the station to the center of the survey network (the average value of the coordinates of all stations), the average division is carried out (standby mode);

Center uniform partition method: divide evenly according to the distance from the site to the center of the network (recommended mode, which can avoid the situation that the baseline length of each area is too large).

In the third step, based on the GNSS data-related products and the station partition list, the parallel automatic GNSS baseline solution for each partition is performed to obtain the single-day baseline solution for each partition.

Daily single-day tasks are independent of each other, and single-day data needs to be partitioned and processed in parallel.

The advantage of the partition + parallel technology lies in: under the premise of ensuring the accuracy, the calculation efficiency is greatly improved, and at the same time, it can meet the one-time fusion calculation of the ultra-large-scale station network.

A separate baseline calculation process is performed on the station data in each sub-division.

For each zone, the principle of GNSS baseline calculation is as follows:

Assuming that satellite i is observed at station j at time epoch _tj , the linearized dual-frequency carrier phase observation equation can be written as:

is the L1 carrier phase observation quantity, the unit is cycle;

is the L2 carrier phase observation quantity, the unit is cycle;

f ₁ is L1 carrier frequency;

f ₂ is the L2 carrier frequency;

is the phase change caused by satellite and receiver clock differences;

k _ij is the effect of ionospheric refraction;

v _ij is the measurement error and the residual error (such as multipath effect) that is not modeled;

where n _ij is the integer ambiguity,

is the initial phase deviation of the receiver,

is the initial phase deviation of the satellite.

The geometric propagation delay is written as:

here

are the three-dimensional geocentric position vectors of the satellite and the station, respectively, and c is the speed of light. The parameters we need to solve are the geocentric position vectors of the stations and satellites.

tropospheric refraction delay time

It is usually calculated from the surface meteorological observation data (or water vapor radiometer WVR data) through the model, and added to the right side of equations (1) and (2) as a phase observation correction. For the residual tropospheric refraction delay part, usually by some zenith Delay parameters (zenith delay parameters) to model and solve together with these parameters when calculating.

We know that the double-difference combination of the phase (dual-frequency phase difference combination forms a new combination) eliminates the phase change caused by the clock difference between the satellite and the receiver, and eliminates the initial phase deviation between the satellite and the receiver, so that the ambiguity has an integer characteristic. Omitting these parts eliminated in the double-difference combination, the phase observation equation is simply written as:

In the formula, n _1ij is the integer number of L1, and n _2ij is the integer number of L2.

If certain constraints are imposed on the effect of ionospheric refraction, an equation can also be listed:

l _ij (t _j )＝k _ij (t _j )/f ₁ +v _kij (t _j ) (6)

In the formula, l _ij (t _j ) is the influence of ionospheric refraction at time t _j . v _kij (t _j ) is v _1ij (t _j ) or v _2ij (t _j ).

In this way, in each observation epoch, the observation equation of the above-mentioned observations is written as an adjustment equation in matrix form, which is convenient for later derivation, calculation and writing; the adjustment equation is:

In the formula,

is the true value of the observation,

is the coefficient matrix,

is the parameter to be estimated,

Correction number for the observational error.

In the above formula are:

A _k =A _k1 =1/f ₁ ; A _k2 =1/f ₂ ; A _k1 =gA _k2 =A _k

The residual vector of double-differenced observations satisfies:

In the above formulas, D is a double-difference operator, which maps all phase one-way observations in the same epoch into independent double-difference observations, A is the linearized coefficient matrix, C is the cofactor matrix, and σ ² is the first The variance of the test unit weight, where the observation weights of L1 and L2 are assumed to be equal. x _a is the unknown parameter to be solved, which includes the coordinate difference of the observation station, the satellite orbit parameter, the pole shift parameter, the ambiguity parameter and the tropospheric zenith delay parameter, etc. x _k is the ionospheric delay parameter.

The fourth step is to check the quality of the baseline solution of each sub-region to obtain a qualified baseline result. If there is any unqualified situation, it is usually necessary to check the quality of the original observation data. If the station is excluded during the calculation, the result of the later stage will not have the result of this station.

The evaluation index of the quality of the baseline results is usually expressed by the normalized root mean square error (Normalized Root Mean Square, NRMS). The standardized root mean square error in the baseline results is used to indicate the degree to which the baseline value calculated by a single period of time deviates from its weighted average value, which is the residual error obtained from the ambiguity solution of the epoch. NRMS is an important index to measure the results of GAMIT solution, the calculation formula is as follows:

In the formula, Y _i is the calculated value of each baseline vector epoch (the coordinate difference between two stations), Y is the true value of the baseline vector, and N is the total number of epochs (that is, how many days of data are calculated in total, and each day is One epoch, each epoch has a solution), n is the number of baselines,

Variance of values solved for each epoch.

The NRMS value should generally be less than 0.3. The smaller the NRMS value, the higher the estimation accuracy of the baseline; otherwise, the lower the accuracy. If the NRMS value is too large (for example, greater than 0.5), it means that the cycle slip may not be completely repaired during the calculation process, or other problems may occur, and further investigation is required. Generally, it is caused by the short observation time of the observation data. It is recommended to increase the observation time of GNSS data (generally not less than 10 hours).

The fifth step is to perform adjustment processing on the single-day baseline solution results of each sub-region obtained from the calculation, and complete the comprehensive adjustment of the subnetwork baseline solution files to obtain the final positioning results of GNSS stations (coordinates of observation stations) .

The principle of the comprehensive adjustment of the baseline solution is as follows:

The output baseline result file provides the coordinates of the stations participating in the calculation, the satellite orbit parameters and the prior value X′ ₀ of the pole shift parameters (other output values of the baseline solution), the prior constraints P _ac of these parameters, and the solution vector δX _c And the covariance matrix C _x of these parameters, which provides favorable conditions for combining single-day solutions using Kalman filtering. In general, these prior constraints are very loose (100 meters for the satellite, 10 meters for the station), and will not affect the result of the solution. If the constraints are so tight that the added value of the ^χ2 test is too large, consider removing the constraints. From the Kalman filter estimation formula, it can be seen that the key issue of Kalman filter estimation is how to construct A _k , y _k , S _k and W _k . We only care about the position and velocity parameters of the station. Therefore, when constructing the above quantities, we ignore the influence of satellite orbit parameters and pole shift parameters, which is beneficial to simplify the formula.

The position and velocity vectors of the station at time t _k are X _k ,

have:

In the above formula, I is the unit matrix. Obviously, the formula (12) is considered as the state transition equation of the station position and velocity parameters, then the state transition matrix is:

If we want to accurately solve the station coordinates and velocity vectors at time t ₀ , let the prior values of coordinates and velocity at this time be X′ ₀ and

Then there are:

y _k is the parameter solution at time t _k , and A _k is the coefficient matrix at time t _k .

The error equation is written as:

δx ₀ ,

_Vk is the correction value of station position and velocity, and can be obtained from the output result file. The next step is to construct the covariance matrix W _k of the random disturbance.

The matrix W _k is a diagonal matrix, and we only need to give the values corresponding to the unknown parameters. For non-random parameters, the value at its corresponding position is zero, and for random parameters, the value at its corresponding position is not zero. Because it can be considered that the drift of some station coordinates obeys the Markov process and is a random walk, so we only need to give its power spectral density PSD value. The corresponding values are:

σ _i ² =PSD×(t _k+1 -t _k ) (17)

is the value of the i-th parameter in W _k . For station coordinates, we define the unit of power spectral density as m ² /year, and a PSD value of 365 is equivalent to a maximum drift of 1 meter per day for the station coordinates. It is very difficult to determine the size of the power spectral density, which can only rely on experience. At the same time, it is also necessary to carefully consider whether the drift of a certain station obeys the Markov process. Only the points with large jumps found through the repeatability test can be treated as random parameters.

The sixth step is to check the quality of the adjustment results obtained from the calculation, so as to obtain the final qualified coordinate results of the GNSS station.

The inspection content mainly judges the errors in the coordinates of the adjustment results. Generally, the errors in the coordinates in the three-dimensional direction are mm level, and the accuracy of the adjustment results is considered to be qualified. If the error exceeds 2cm, it is judged that the accuracy of the adjustment result is unqualified. In case of non-conformity, it is usually necessary to check the data quality of the original observation data.

In the example provided by the present invention, the large-scale GNSS network automatic parallel solution method based on dynamic partitioning automatically downloads the IGS stations around the GNSS network and the products of the International Analysis Center, and passes the measuring stations covered by the GNSS network through a certain A partitioning method (geographical area partitioning method, central diffusion partitioning method, central uniform partitioning method) is automatically divided into several sub-partitions, and then based on multi-core parallel technology, parallel baseline calculations are performed on each partition, and then based on the baseline solution equation file of each partition Carry out sub-network combined adjustment, and GNSS positioning results can be obtained through comprehensive adjustment. Finally, the baseline solution and adjustment solution can be evaluated and checked for accuracy. It does not require user interaction, and can automatically complete a series of tasks such as data download, subnet partitioning, baseline calculation, subnet fusion, and comprehensive adjustment, thereby reducing server pressure, and the algorithm is not limited by the number of GNSS stations. Parallel computing technology greatly improves the efficiency of GNSS data computing, and provides an algorithm basis for comprehensive positioning services of large-scale GNSS networks.

Example 2

Based on the same inventive concept as in Embodiment 1, a large-scale GNSS network parallel computing system based on dynamic partitioning of the present invention includes:

A data acquisition module for collecting observation data;

The partition processing module is used to partition the observation stations corresponding to the observation data into subnetworks to obtain a list of station partitions;

For the specific implementation scheme of each module in the system of the present invention, refer to the specific implementation process of each step in the method.

Those skilled in the art should understand that the embodiments of the present application may be provided as methods, systems, or computer program products. Accordingly, the present application may take the form of an entirely hardware embodiment, an entirely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, the present application may take the form of a computer program product embodied on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) having computer-usable program code embodied therein.

The present application is described with reference to flowcharts and/or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of the present application. It should be understood that each procedure and/or block in the flowchart and/or block diagram, and a combination of procedures and/or blocks in the flowchart and/or block diagram can be realized by computer program instructions. These computer program instructions may be provided to a general purpose computer, special purpose computer, embedded processor, or processor of other programmable data processing equipment to produce a machine such that the instructions executed by the processor of the computer or other programmable data processing equipment produce a An apparatus for realizing the functions specified in one or more procedures of the flowchart and/or one or more blocks of the block diagram.

These computer program instructions may also be stored in a computer-readable memory capable of directing a computer or other programmable data processing apparatus to operate in a specific manner, such that the instructions stored in the computer-readable memory produce an article of manufacture comprising instruction means, the instructions The device realizes the function specified in one or more procedures of the flowchart and/or one or more blocks of the block diagram.

These computer program instructions can also be loaded onto a computer or other programmable data processing device, causing a series of operational steps to be performed on the computer or other programmable device to produce a computer-implemented process, thereby The instructions provide steps for implementing the functions specified in the flow diagram procedure or procedures and/or block diagram procedures or blocks.

The above is only a preferred embodiment of the present invention, it should be pointed out that for those of ordinary skill in the art, without departing from the technical principle of the present invention, some improvements and modifications can also be made, these improvements and modifications It should also be regarded as the protection scope of the present invention.

Claims

A large-scale GNSS network parallel solution method based on dynamic partitioning, characterized in that it includes:

collect observational data;

Carry out subnetwork partitioning for the observation stations corresponding to the observation data to obtain a list of station partitions;

Based on the station partition list and observation data, parallel baseline calculations are performed for each partition to obtain a single-day baseline solution for each partition;

Adjustment processing is carried out on the single-day baseline solution of each division to obtain the final positioning results of each observation station.
A kind of large-scale GNSS network parallel solution method based on dynamic partition according to claim 1, it is characterized in that, after described collecting observation data, also comprise

The observation data is preprocessed to obtain complete and qualified observation data products. The preprocessing includes data continuity check, integrity check, data editing and data quality check processing.
A kind of large-scale GNSS network parallel solution method based on dynamic partitioning according to claim 1, it is characterized in that, described subnetwork partitioning method comprises geographic region partitioning method, central diffusion partitioning method and center uniform partitioning method any one kind.
A kind of large-scale GNSS network parallel solution method based on dynamic partition according to claim 1, it is characterized in that, described baseline solution is carried out to each partition, obtains each partition single-day baseline solution, comprising:

Assuming that satellite i is observed at station j at time epoch tj , the linearized dual-frequency carrier phase observation equation is written as:

Among them, the subscripts 1 and 2 represent the carrier phases of L1 and L2 respectively;

is the L1 carrier phase observation quantity, the unit is cycle;

is the L2 carrier phase observation quantity, the unit is cycle;

f 1 is L1 carrier frequency;

f 2 is the L2 carrier frequency;

τ ij is the geometric propagation delay time of the carrier signal between the satellite and the receiver;

is the phase change caused by satellite and receiver clock differences;

is the tropospheric refraction delay on the propagation path of the carrier signal;

k ij is the effect of ionospheric refraction;

v ij is measurement error and residual error;

where n ij is the integer ambiguity,
is the initial phase deviation of the receiver,
is the initial phase deviation of the satellite;

The geometric propagation delay is written as:

here
are the three-dimensional geocentric position vectors of the satellite and the station respectively, and c is the speed of light; the parameters to be solved are the geocentric position vectors of the station and the satellite;

If the phase change eliminated in the double-difference combination is omitted, the phase observation equation is simply written as:

In the formula, n 1ij is the integer number of L1, n 2ij is the integer number of L2;

If certain constraints are imposed on the effect of ionospheric refraction, an equation is also listed:

l ij (t j )＝k ij (t j )/f 1 +v kij (t j ) (6)

In the formula, l ij (t j ) is the influence of ionospheric refraction at time t j ; v kij (t j ) is v 1ij (t j ) or v 2ij (t j );

In each observation epoch, the observation equation of the above observations is written as the adjustment equation in matrix form:

In the formula,
is the true value of the observation,
is the coefficient matrix,
is the parameter to be estimated,
is the observation error correction number;

In the above formula are:

A k =A k1 =1/f 1 ; A k2 =1/f 2 ; A k1 =gA k2 =A k

The residual vector of double-differenced observations satisfies:

In the formula, D is a double-difference operator, which maps all phase one-way observations in the same epoch into independent double-difference observations, A is the linearized coefficient matrix, C is the cofactor matrix, and σ 2 is the prior unit weight Variance, x a is the unknown parameters to be solved, including the coordinate difference of observation stations, satellite orbit parameters, pole shift parameters, ambiguity parameters and tropospheric zenith delay parameters, x k is the ionospheric delay parameters.
A kind of large-scale GNSS network parallel solution method based on dynamic partition according to claim 1, it is characterized in that, after described obtaining each partition single-day baseline solution, also include:

Perform a quality check on the baseline solution of each partition to obtain a qualified baseline result;

If there are unqualified conditions, check the quality of the original observation data. If it is caused by the poor quality of the station data, remove the station during the calculation.
A kind of large-scale GNSS network parallel solution method based on dynamic partition according to claim 5, it is characterized in that, the concrete process of described quality inspection is:

Calculate the evaluation index of the baseline result, the evaluation index is expressed by the standardized root mean square error NRMS, and the calculation formula is as follows:

In the formula, Y i is the calculated value of each baseline vector epoch, Y is the true value of the baseline vector, N is the total number of epochs, n is the number of baselines,
is the variance of the calculated value for each epoch;

If the NRMS value is less than the set threshold, it is judged that the baseline solution is qualified, otherwise it is not qualified.
A kind of large-scale GNSS network parallel solution method based on dynamic partition according to claim 1, it is characterized in that, described each partition single-day baseline solution is carried out adjustment processing, obtains the final positioning result of each observation station, include:

Let the position and velocity vectors of the observation station at time t 0 be X 0 ,
The position and velocity vectors of the station at time t k are X k ,
Then there are:

The third item on the right side of the above formula is the random drift part of the station between t 0 and t k , r is the random drift coefficient, δρ k is the observation error at each moment, these influences are ignored, and the speed of the station is assumed does not change with time, formula (11) is written in matrix form as follows:

In the above formula, I is the unit matrix, and the formula (12) is considered as the state transition equation of the station position and velocity parameters, then the state transition matrix is:

Let the prior values of coordinates and velocity at this moment be X′ 0 and
Then there are:

y k is the parameter solution at time t k , and A k is the coefficient matrix at time t k ;

The error equation is written as:

δx 0 ,
is the correction amount of station position and velocity, V k is obtained from the output result file;

The covariance matrix W k of the random disturbance is a diagonal matrix. Since the drift of the station coordinates is subject to the Markov process, it is only necessary to give its power spectral density PSD value; the corresponding value is:

σ i 2 =PSD×(t k+1 -t k ) (17)

is the value of the i-th parameter in W k , PSD is the empirical value;

Based on the constructed A k , y k , S k and W k , the Kalman filter estimation algorithm is used for adjustment processing, and the final positioning results of each observation station are obtained.
A kind of large-scale GNSS network parallel solution method based on dynamic partition according to claim 1, it is characterized in that, after the final positioning result of each observation station is obtained, it also includes:

Perform quality checks on the final positioning results of each observation station to obtain the final qualified coordinate results of the GNSS station.
A large-scale GNSS network parallel computing system based on dynamic partitioning, characterized in that it includes:

A data acquisition module for collecting observation data;

The partition processing module is used to partition the observation stations corresponding to the observation data into subnetworks to obtain a list of station partitions;

The parallel calculation module is used to perform parallel baseline calculations for each zone based on the station zone list and observation data, and obtain a single-day baseline solution for each zone;

And, the comprehensive adjustment module is used for adjusting the single-day baseline calculation of each subregion to obtain the final positioning results of each observation station.