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

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

Info

Publication number
WO2022048694A1
WO2022048694A1 PCT/CN2021/123845 CN2021123845W WO2022048694A1 WO 2022048694 A1 WO2022048694 A1 WO 2022048694A1 CN 2021123845 W CN2021123845 W CN 2021123845W WO 2022048694 A1 WO2022048694 A1 WO 2022048694A1
Authority
WO
WIPO (PCT)
Prior art keywords
formula
value
iteration
station
satellite
Prior art date
Application number
PCT/CN2021/123845
Other languages
English (en)
French (fr)
Inventor
郭金运
郭恒洋
杨洲铭
邢云鹏
刘新
孔巧丽
Original Assignee
山东科技大学
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 山东科技大学 filed Critical 山东科技大学
Publication of WO2022048694A1 publication Critical patent/WO2022048694A1/zh
Priority to ZA2022/03722A priority Critical patent/ZA202203722B/en
Priority to US17/730,567 priority patent/US20220299652A1/en

Links

Images

Classifications

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

Definitions

  • the invention belongs to the technical field of satellite navigation and positioning, and relates to a GNSS single-point positioning method based on spherical harmonic expansion.
  • the GNSS single-point positioning method refers to the position of the satellite in space and the satellite clock difference at the moment of observation given by the satellite ephemeris, and the distance from the satellite to the GNSS receiver measured by a GNSS receiver.
  • the rendezvous method is a positioning method for independently determining the three-dimensional coordinates of the GNSS receiver in the earth coordinate system.
  • the positioning method using the satellite position and satellite clock offset provided by the broadcast ephemeris, as well as the pseudo-range observations, is called standard single-point positioning. Due to the limitation of the accuracy of broadcast ephemeris and pseudorange observations, and the influence of various correction errors in the mathematical model cannot be completely eliminated, the accuracy of standard single-point positioning can usually only reach the meter level. Therefore, standard single-point positioning is mainly used in low-precision positioning fields such as navigation positioning, resource survey, and survey. Using precise ephemeris and precise satellite clock error data, the positioning method based on phase observation data is called precise single-point positioning, and the precision of precise single-point positioning can reach the order of centimeters.
  • the errors related to signal propagation (mainly referring to ionospheric delay and tropospheric delay errors) cannot be well modeled, and therefore have a greater impact on positioning.
  • the single-frequency receiver needs to be corrected by the model; for the dual-frequency receiver, the influence of the ionospheric delay can be eliminated by the linear combination method; for the tropospheric delay, due to the non-dispersive nature of the troposphere, the dual-frequency receiver cannot be used.
  • the method of linear combination of observation data eliminates tropospheric delay error, and its influence on positioning accuracy cannot be ignored.
  • the tropospheric delay is divided into two parts: dry delay and wet delay.
  • dry delay in the zenith direction is calculated by the tropospheric delay model, and it is mapped to the propagation path through the mapping function; while the wet delay in the zenith direction is solved as an unknown parameter, and through the mapping function map to the propagation path.
  • the tropospheric delay models or methods in GNSS measurements include the Sastamonen model, the UNB3 model, and the meteorological element method.
  • the above three tropospheric delay correction methods have the following problems: the accuracy of the Sastamonen model is relatively high, but it needs to use the measured meteorological data, and the acquisition of meteorological data limits the application of the Sastamonen model; UNB3 The model does not require the injection of measured meteorological data, but the accuracy of the UNB3 model is limited; the meteorological element method needs to use external meteorological equipment to collect data, the accuracy of different equipment is different, and the equipment is cumbersome and expensive, which increases the workload of field data collection; The meteorological element method also reduces the work efficiency of single-point positioning, and it is difficult to provide the required data for all-weather single-point positioning.
  • tropospheric delay correction methods build a variety of models to meet the needs of tropospheric delay correction under different conditions. Considering the practicability and universality, the workload of the measurement field and office will increase, and the data processing efficiency will decrease.
  • the tropospheric delay error correction mostly uses empirical models, and the parameter estimation of such empirical models also mostly uses empirical values, so it is difficult to conform to the parameter estimation in actual single-point positioning.
  • the purpose of the present invention is to propose a GNSS standard single-point positioning method based on spherical harmonic expansion, so that the standard single-point positioning can be performed quickly, efficiently and accurately.
  • a GNSS standard single-point positioning method based on spherical harmonic expansion is used to describe the error terms related to the altitude and azimuth between the station and the satellite, including the tropospheric delay error and the ionospheric delay error, using GNSS observations.
  • the data and satellite ephemeris data can realize the calculation of the spatial position of the station, and then realize the GNSS standard single-point positioning.
  • the method of the invention can obtain the position information of the measuring point efficiently and quickly, and has the advantages of simple operation, convenient data processing and high calculation efficiency.
  • Embodiment 1 is a schematic flowchart of the GNSS standard single-point positioning based on spherical harmonic expansion in Embodiment 1 of the present invention
  • FIG. 2 is a schematic flow chart of GNSS precise single-point positioning based on spherical harmonic expansion in Embodiment 2 of the present invention
  • FIG. 3 is a statistical diagram of the results of a GNSS standard single-point positioning method based on spherical harmonic expansion according to an embodiment of the present invention
  • FIG. 4 is a statistical diagram of the results of a GNSS precise single-point positioning method based on spherical harmonic expansion according to an embodiment of the present invention.
  • This embodiment describes a GNSS standard single-point positioning method based on spherical harmonic expansion.
  • the GNSS standard single-point positioning method based on spherical harmonic expansion includes the following steps:
  • the satellite position corresponding to each epoch is obtained by interpolation of the broadcast ephemeris; and the satellite clock error is obtained by using the satellite clock error parameters provided by the broadcast ephemeris.
  • ground receiver clocks use quartz clocks.
  • the accuracy of ground receivers is lower than that of satellite clocks, and the clock error changes irregularly and cannot be well fitted by polynomials. Therefore, the receiver clock error tr is used as the parameter to be estimated for estimation.
  • i represents the i-th observation epoch
  • j represents the satellite number
  • ⁇ represents the pseudorange observation value
  • ⁇ represents the pseudorange observation value of satellite j at the i-th observation epoch
  • Re represents the geometric distance between the position of the receiver and the position of satellite j at the i-th observation epoch
  • c is the speed of light
  • tr is the station receiver clock error
  • ( t r ) i is the i-th observation epoch.
  • Receiver clock error represents the satellite clock error, represents the clock error of satellite j at the i-th observation epoch; represents the tropospheric delay error of satellite j on the signal propagation path of the i-th observation epoch; represents the ionospheric delay error of satellite j on the signal propagation path of the i-th observation epoch;
  • ⁇ ⁇ represents the residual error of the pseudo-range observation data;
  • the spatial three-dimensional coordinates of the satellite j at the i-th observation epoch at the moment of observation are defined as
  • the approximate coordinates of the station are (X 0 , Y 0 , Z 0 ), then the geometric distance from the satellite to the approximate position of the station Expressed as:
  • GNSS pseudorange observation data types are divided into single-frequency observation data and dual-frequency observation data. If the GNSS receiver for standard single-point positioning is a dual-frequency receiver, the ionospheric delay error is eliminated by the combination of dual-frequency ionosphere cancellation, and the combination is as follows: ⁇ is the combined observation value of pseudorange and ionosphere; f 1 , f 2 are the frequencies of the corresponding observation value; ⁇ 1 , ⁇ 2 respectively represent the pseudorange observation value corresponding to the two frequencies.
  • the first-order main term of the ionospheric delay error is eliminated by the combination of dual-frequency elimination of the ionosphere, and the ionospheric delay error term is 0.
  • the spherical harmonic expansion is mainly used to represent the tropospheric delay error term, and the order of the spherical harmonic expansion is different from the observation equation constructed by the single-frequency pseudorange observations.
  • both tropospheric delay error and ionospheric delay error are a kind of signal delay error generated when the satellite signal passes through the atmosphere, and both are related to the propagation of the satellite signal; and the propagation path of the satellite signal is also It is related to the position of the station and the satellite.
  • the altitude and azimuth between the station and the satellite can be calculated by using the coordinates of the station and the satellite.
  • spherical harmonic expansion the altitude and the satellite can be expressed.
  • Azimuth-related error terms including tropospheric delay error and ionospheric delay error, to eliminate or minimize their impact on standard single-point positioning results.
  • the standard single-point positioning observation equation based on spherical harmonic expansion is expressed as:
  • n is the order of spherical harmonic expansion
  • m is the order of spherical harmonic expansion
  • Nmax is the maximum order of spherical harmonic expansion
  • C nm and S nm respectively represent the coefficients of the spherical harmonic function model of order n of order m
  • C nm and S nm are the parameters to be obtained for spherical harmonic expansion
  • the standard single-point positioning observation equation based on spherical harmonic expansion expressed by formula (2) does not involve the correction model of tropospheric delay error, and does not need to consider the application of different models in different regions, and in this embodiment, the coefficient of spherical harmonic expansion is regarded as
  • the parameters to be found are directly solved to correct the error terms related to the altitude and azimuth angles between the station and the satellite, mainly including the tropospheric delay error, without considering the meteorological information, which is more universal.
  • v ⁇ represents the correction value of the pseudorange observation data
  • the parameters to be determined include 3 position parameters (dX, dY, dZ), epoch receiver clock error dtr parameters and (Nmax+1) 2 spherical harmonic expansion parameters C nm and S nm .
  • the optimal solution of all unknowns can be obtained by least squares. Therefore, as long as the correction term can be described linearly, the number of observation equations can be increased by increasing the observation epoch to solve these unknown parameters.
  • this embodiment selects an appropriate sliding window, and performs sliding calculation on the simplified linear expression (6) of the standard single-point positioning error equation.
  • X [dX dY dZ (dt r ) 1 ... (dt r ) i C 00 ... C NmaxNmax S NmaxNmax ] T .
  • V BX-L (7)
  • P is the unit weight matrix.
  • the empirical threshold of the condition number is set through the calculated condition number of the coefficient matrix B, and the condition number is compared with the empirical threshold. When the condition number is greater than the empirical threshold, it indicates that the error equation (7) is ill-conditioned, otherwise it is Well-being.
  • this embodiment proposes to use the method of truncated singular value decomposition to provide the initial value of the parameter X to be determined for the least squares spectrum correction iteration, and calculate the to be determined based on the least squares spectrum correction iteration The value of parameter X.
  • the idea of the truncated singular value method is to calculate the mean of the singular values in the matrix S, and use the mean as the threshold for truncating singular values; filter out the singular values smaller than the threshold and zero them.
  • the truncated singular value decomposition method converts the least squares method of the equation into a direct solution method, the solution of the ill-conditioned equation is avoided to deviate too much from the real solution, and the ill-conditioned problem of the coefficient matrix B of the error equation is effectively solved.
  • coefficient matrix B R n ⁇ m , R n ⁇ m represents a real number matrix with n rows and m columns; then the singular value decomposition of coefficient matrix B is:
  • the truncated singular value matrix Bk of the coefficient matrix B ⁇ Rn ⁇ m is defined as:
  • K is an arbitrary real number.
  • the least squares spectrum correction iteration is an iterative calculation method based on least squares; iteratively calculates the solution of the unknown parameter X based on the least squares spectrum correction iteration.
  • the solution process is as follows:
  • the iterative solution X based on least-squares spectrum correction includes two iterative processes: iteration 1 and iteration 2; wherein, a first comparison threshold is set for judging whether iteration 1 converges, and a second comparison threshold for judging whether iteration 2 converges is set Compare thresholds.
  • Iteration 1 In the first iteration, set the initial value of K to 1, and solve the truncated singular value obtained by formula (11). As the initial value of X in the first iteration of formula (14), and assign it to X on the right side of the equation of formula (14);
  • the difference obtained after the difference operation is greater than the first comparison threshold, it means that the iteration does not converge, and the K value in the least squares spectrum correction iteration needs to be corrected, that is, the K value is incremented by 1, and the above iteration 1 process is continued;
  • the iteration 1 is jumped out and the process of iteration 2 is entered.
  • Iteration 2 compares the value of the unknown parameter X obtained by using the formula (14) with the second comparison threshold;
  • the position parameters (dX, dY, dZ) are added to the approximate coordinates of the station (X 0 , Y 0 , Z 0 ), update the approximate coordinates of the station, linearize the standard single-point positioning error equation after updating the approximate coordinates of the station by formula (5), and simplify the linearized equation by formula (6).
  • the value of the unknown parameter X obtained by using the formula (14) is less than or equal to the second comparison threshold, it means that the iteration 2 converges and the iteration is jumped out; at this time, the value of the unknown parameter X obtained by using the formula (14) is the unknown parameter X. solution.
  • step I.3.4 the process of solving the solution of the unknown parameter X based on the least squares method is as follows:
  • the position parameters (dX, dY, dZ) are added to the approximate coordinates of the station (X 0 , Y 0 , Z 0 ), update the approximate coordinates of the station, perform the calculation process of formula (5)-formula (8), and re-solve the value of the unknown parameter X;
  • the value of the unknown parameter X obtained by using the formula (8) is less than or equal to the third comparison threshold, it means that the iteration has converged and the iteration is skipped; at this time, the value of the unknown parameter X obtained by using the formula (8) is the unknown parameter. solution for X.
  • the parameter values contained in the unknown parameter X are obtained, that is, the position parameters (dX, dY, dZ) of the station, the receiver clock error dt r and the coefficient of spherical harmonic expansion C nm and S nm ;
  • the position parameters (dX, dY, dZ) in the unknown parameter X are respectively added to the approximate coordinates (X 0 , Y 0 , Z 0 ) of the station, that is, in the geocentric geofixed coordinate system obtained by standard single-point positioning The station coordinates (X, Y, Z) of .
  • the specific application process is: reading the existing spherical harmonic expansion coefficients, and expressing the error terms related to the altitude and azimuth angles between the station and the satellite through spherical harmonic expansion, mainly including tropospheric delay errors and ionospheric delay errors;
  • spherical harmonic expansion mainly including tropospheric delay errors and ionospheric delay errors
  • the single-point positioning observation equation is used, the correction of the error terms related to the altitude and azimuth angles between the station and the satellite is no longer considered, and the solving process of the tropospheric delay error and ionospheric delay error in positioning is omitted, so that the positioning equation is only It is only necessary to calculate the position parameters of the station and the clock error parameters of the station receiver, which greatly simplifies the standard single-point positioning process.
  • I.3.6 Convert the coordinates of the measuring station from the earth-centered ground-fixed coordinate system to the station-centered coordinate system to realize GNSS standard single-point positioning.
  • (X, Y, Z) are the coordinates of the fixed coordinate system at the center of the earth, (lon, lat, alt) are the coordinates of the latitude and longitude coordinate system; M is the radius of curvature of the reference ellipsoid, and E is the eccentricity of the ellipsoid.
  • the corresponding coordinates of the station’s approximate coordinates (X 0 , Y 0 , Z 0 ) in the latitude, longitude and altitude coordinate system are calculated by formula (15), That is (lon 0 ,lat 0 ,alt 0 ).
  • (e 1 , n 1 , u 1 ) is the position of (X, Y, Z) coordinates obtained by taking the approximate coordinates of the station (X 0 , Y 0 , Z 0 ) as the origin.
  • the station coordinates are converted from the geocentric ground-fixed coordinate system to the station center coordinate system.
  • the GNSS standard single-point positioning method in this embodiment 1 represents the error terms related to the azimuth and elevation angles between the station and the satellite based on spherical harmonic expansion, so that standard single-point positioning can be performed quickly, efficiently and accurately.
  • This embodiment 2 describes a GNSS precise single-point positioning method based on spherical harmonic expansion.
  • the GNSS precise single-point positioning method based on spherical harmonic expansion includes the following steps:
  • the observation data broadcast ephemeris data and precise ephemeris data, preprocess the carrier phase observation data, and calculate the errors that are irrelevant or not closely related to the altitude and azimuth between the station and the satellite; the above errors mainly include calculation Satellite clock error, detection cycle slip and initial value of carrier phase integer cycle ambiguity parameters are calculated.
  • the accuracy of broadcast ephemeris cannot meet the requirements of precise single-point positioning, so the satellite precise ephemeris must be used to calculate the satellite position; the satellite clock error estimated by second-order polynomial cannot meet the accuracy requirements of precise single-point positioning; Since the clock errors of different satellites are different, their magnitude must be determined in advance, and then substituted into the observation equation to eliminate its influence.
  • the post-event precise ephemeris with sampling interval of 15min and 5min and satellite clock error products with sampling interval of 5min and 30s provided by IGS can meet the precision requirements of precise positioning.
  • the sampling interval of the observed values is generally smaller than the sampling interval of the data in the above files. Therefore, it is necessary to interpolate and calculate the satellite position and satellite clock error corresponding to each epoch.
  • i represents the i-th observation epoch
  • j represents the satellite number
  • tr represents the receiver clock error
  • ( t r ) i represents the receiver clock error at the i-th observation epoch
  • ts represents the satellite clock Difference, represents the clock error of satellite j at the i-th observation epoch
  • GNSS carrier phase observation data types are divided into single-frequency observation data and dual-frequency observation data.
  • the ionospheric delay error can be eliminated by the combination of dual-frequency ionosphere cancellation.
  • the specific combination is as follows: In the formula, represents the combined observation value of carrier phase cancellation ionosphere (equivalent distance); f 1 , f 2 represent the frequency of the corresponding observation value; respectively represent the carrier phase observations (equivalent distances) corresponding to the two frequencies. If the observation data is a dual-frequency observation value, in the above formula (1), the first-order main term of the ionospheric delay error is eliminated by the combination of dual-frequency elimination of ionosphere, and the ionospheric delay error term is 0.
  • the spherical harmonic expansion is mainly used to represent the tropospheric delay error term, and the order of the spherical harmonic expansion is different from the observation equation constructed by the single-frequency carrier phase observations.
  • the carrier phase observations are used to form GF (geometry-independent) combined observations to solve the ambiguity of the whole cycle and detect the cycle slip.
  • GF geometry-independent
  • the formula represents the observed value of the GF combined carrier phase (equivalent distance); respectively represent the observed values of the carrier phase corresponding to the two frequencies; ⁇ 1 and ⁇ 2 represent the wavelengths of the carrier phases corresponding to the two frequencies respectively; f 1 , f 2 represent the frequencies of the two carrier phases.
  • the GF combined observations are independent of geometric distance, so this combination is suitable for detecting cycle slips.
  • the error terms related to the altitude and azimuth angles between the station and the satellite are expressed, wherein the error terms mainly include the tropospheric delay error and the ionospheric delay error.
  • the tropospheric delay error and ionospheric delay error are both a kind of signal delay error generated when the satellite signal passes through the atmosphere, and both are related to the propagation of the satellite signal; and the propagation path of the satellite signal is related to the station and The position of the satellite is related.
  • the altitude and azimuth between the station and the satellite can be calculated by using the coordinates of the station and the satellite.
  • the spherical harmonic expansion can be used to express the altitude and azimuth between the station and the satellite.
  • the error term mainly includes tropospheric delay error and ionospheric delay error, in order to eliminate or minimize its influence on precise single-point positioning results.
  • the precise single-point positioning observation equation after spherical harmonic expansion is expressed as:
  • n is the order of spherical harmonic expansion
  • m is the order of spherical harmonic expansion
  • Nmax is the maximum order of spherical harmonic expansion
  • P nm (cose) represents the associative Legendre polynomial of order n of order m
  • C nm and S nm represent the coefficients of spherical harmonic expansion of n order m order
  • the precise single-point positioning observation equation based on the spherical harmonic expansion proposed by formula (2) does not involve the correction model of the tropospheric delay error, and does not need to consider the application of different models in different regions.
  • the coefficient of the spherical harmonic expansion is taken as The parameters to be found are directly solved to correct the error terms related to the altitude and azimuth angles between the station and the satellite, mainly including the tropospheric delay error, without considering the meteorological information, which is more universal; for the convenience of expressing the spherical harmonic function expansion ,make:
  • formula (5) represents the initial ambiguity value of satellite j at the i-th epoch; represents the integer ambiguity parameter of the phase observation data of satellite j at the i-th epoch, dN represents the correction number of the integer ambiguity parameter of the phase observation data, Represents the correction number of the integer ambiguity parameter of the phase observation data of satellite j at the i-th epoch.
  • the correction number is used as the number to be determined.
  • the parameters are solved, and when solving as the parameters to be solved,
  • the coefficient of is 1;
  • d(t r ) i represents the parameter to be determined of the receiver clock error of the station at the i-th observation epoch, which is similar to the processing method of the integer ambiguity, because the receiver clock error of the station cannot be well It is fitted by a polynomial, so it is solved as a parameter to be solved.
  • the coefficient of dt r is 1.
  • the parameters to be determined include: 3 position parameters (dX, dY, dZ), epoch receiver clock error parameters dt r , j ambiguity parameters dN and (Nmax+1) 2 spherical harmonic expansion coefficients C nm , S nm .
  • the optimal solution of all unknowns can be obtained by least squares.
  • the number of observation equations can be increased by increasing the observation epoch to solve these unknown parameters.
  • V BX-L (7)
  • the GNSS precise single-point positioning method based on spherical harmonic expansion is based on this idea, and the spherical harmonic expansion is used to represent the altitude and azimuth between the station and the satellite.
  • the error term of directly solves the coefficients of spherical harmonic expansion in the precise single-point positioning equation.
  • P is the unit weight matrix.
  • condition number of the coefficient matrix B of the precise single-point positioning error equation is calculated, the empirical threshold of the condition number is set, and the condition number is compared with the empirical threshold. When the condition number is greater than the empirical threshold, the error equation is ill-conditioned.
  • this embodiment proposes a method of using truncated singular value decomposition to provide the initial value of the parameter X to be determined for the least squares spectrum correction iteration, and calculate the to-be-required parameter X based on the least squares spectrum correction iteration. Find the value of parameter X. Since the truncated singular value decomposition method converts the least squares method of the equation into a direct solution method, the solution of the ill-conditioned equation is avoided to deviate too much from the real solution, and the ill-conditioned problem of the coefficient matrix B of the error equation is effectively solved.
  • the truncated singular value matrix Bk of the coefficient matrix B ⁇ Rn ⁇ m is defined as:
  • the least squares spectral correction iteration is an iterative calculation method based on least squares.
  • K is an arbitrary real number; the solution of the unknown parameter X is calculated iteratively based on the least squares spectral correction.
  • the iterative solution X based on least-squares spectrum correction includes two iterative processes: iteration 1 and iteration 2; wherein, a first comparison threshold is set for judging whether iteration 1 converges, and a second comparison threshold for judging whether iteration 2 converges is set comparison threshold;
  • Iteration 1 In the first iteration, set the initial value of K to 1, and solve the truncated singular value obtained by formula (11). As the initial value of X in the first iteration of formula (14), and assign it to X on the right side of the equation of formula (14);
  • the difference obtained after the difference operation is greater than the first comparison threshold, it means that the iteration does not converge, and the K value in the least squares spectrum correction iteration needs to be corrected, that is, the K value is incremented by 1, and the above iteration 1 process is continued;
  • Iteration 2 compares the value of the unknown parameter X obtained by using the formula (14) with the second comparison threshold;
  • the position parameters (dX, dY, dZ) are added to the approximate coordinates of the station (X 0 , Y 0 , Z 0 ), update the approximate coordinates of the station, linearize the precision single-point positioning error equation after updating the approximate coordinates of the station by formula (5), and simplify the precision after linearization by formula (6).
  • the value of the unknown parameter X obtained by using the formula (14) is less than or equal to the second comparison threshold, it means that the iteration 2 has converged, and the iteration is jumped out; at this time, the value of the unknown parameter X obtained by using the formula (14) is the value to be solved. Solution for unknown parameter X.
  • step II.3.4 Use the least squares method to solve for the unknown parameter X, go to step II.3.5.
  • step II.3.4 the process of solving the solution of the unknown parameter X based on the least squares method is as follows:
  • the position parameters (dX, dY, dZ) are added to the approximate coordinates of the station (X 0 , Y 0 , Z 0 ), update the approximate coordinates of the station;
  • the value of the unknown parameter X obtained by using the formula (8) is less than or equal to the third comparison threshold, it means that the iteration has converged and the iteration is skipped; at this time, the value of the unknown parameter X obtained by using the formula (8) is the solution to be solved.
  • the position parameters (dX, dY, dZ) are respectively added to the approximate coordinates (X 0 , Y 0 , Z 0 ) of the station, that is, the station coordinates ( X, Y, Z).
  • the obtained spherical harmonic expansion coefficients C nm and S nm can be directly used for subsequent precise single-point positioning.
  • the specific application process is as follows:
  • the GNSS precise single-point positioning method described in the second embodiment is based on spherical harmonic expansion to express the error terms related to the azimuth and altitude between the station and the satellite, so that precise single-point positioning can be performed quickly, efficiently and accurately .
  • the sampling interval of GNSS observation data is 30s.
  • the observation data is selected from dual-frequency pseudorange observations and dual-frequency carrier phase observations of GPS satellites.
  • the reference coordinates of the site are the site coordinates provided by IGS.
  • the observation data selects the GPS satellite dual-frequency pseudorange observation value, then the ionospheric delay error is eliminated by the pseudorange linear combination.
  • the error terms related to the altitude and azimuth angles between the satellite and the station are represented by spherical harmonic expansion, mainly the tropospheric delay error.
  • the coefficients of each order of spherical harmonic expansion and the receiver clock error of the station need to be solved as parameters to be determined.
  • the parameters to be determined include: correction values of point coordinates in 3 directions, 10 receiver clock errors (one receiver clock error parameter per epoch), and 16 coefficients of spherical harmonic expansion to the third order: A 00 , A 10 , A 11 , B 11 , A 20 , A 21 , B 21 , A 22 , B 22 , A 30 , A 31 , B 31 , A 32 , B 32 , A 33 , B 33 .
  • the cutoff elevation angle of the satellite is set to 5°, and the coefficient k of the least squares spectral correction iteration is set to 2.
  • the RMS in the E direction of the positioning result calculated by the GNSS high-precision data processing software Bernese software is better than 1.28m
  • the RMS in the N direction is better than 1.61m
  • the RMS in the U direction is better than 3.14m.
  • the observation data selects the GPS satellite dual-frequency carrier phase observation value, and the ionospheric delay error is eliminated by the linear combination of the carrier phase.
  • the error terms related to the altitude and azimuth angles between the satellite and the station are represented by the spherical harmonic expansion, mainly the tropospheric delay error; the coefficients of each order of the spherical harmonic expansion need to be solved as the parameters to be determined.
  • the parameters to be determined include: correction values of point coordinates in 3 directions, 15 receiver clock errors (one receiver clock error parameter per epoch), ambiguity parameters of each satellite in each epoch, spherical harmonic expansion 49 coefficients up to 6th order 6th degree: A 00 , A 10 , A 11 , B 11 , A 20 , A 21 , B 21 , A 22 , B 22 , A 30 , A 31 , B 31 , A 32 , B 32 , A33 , B33 , A40 , A41 , B41 , A42 , B42 , A43 , B43 , A44 , B44 , A50 , A51 , B51 , A52 , B52 , A 53 , B 53 , A 54 , B 54 , A 55 , B 55 , A 60 , A 61 , B 61 , A 62 , B 62 , A 63 , B 63 , A 64 , B 64 , A
  • the cut-off altitude angle of the satellite is set to 5°, and the coefficient k of the least squares spectral correction iteration is set to 3.
  • Table 4 shows the standard single-point positioning results calculated by the internationally renowned GNSS high-precision data processing software Bernese software developed by the Institute of Astronomy, University of Bern, Switzerland. It can be seen that the RMS in the E direction of the positioning result calculated by the GNSS high-precision data processing software Bernese software is better than 0.009m, the RMS in the N direction is better than 0.02m, and the RMS in the U direction is better than 0.03m.
  • the present invention is based on the GNSS single-point positioning method of spherical harmonic expansion, and uses observation data and satellite ephemeris data to efficiently and quickly obtain the point information of the measuring point. High, no need to provide meteorological documents, easy operation, low work difficulty, reliable measurement results, and high measurement accuracy.

Landscapes

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

Abstract

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

Description

一种基于球谐展开的GNSS单点定位方法 技术领域
本发明属于卫星导航定位技术领域，涉及一种基于球谐展开的GNSS单点定位方法。
背景技术
GNSS单点定位方法，是指根据卫星星历给出的观测瞬间卫星在空间的位置和卫星钟差，以及由一台GNSS接收机所测定的从卫星到GNSS接收机之间的距离，通过距离交会的方法，独立测定该GNSS接收机在地球坐标系中的三维坐标的定位方法。
利用广播星历提供的卫星位置和卫星钟差，以及伪距观测值进行的定位方式，称为标准单点定位。由于广播星历和伪距观测值精度的限制，且其数学模式中各项改正误差的影响又无法完全消除，因而，标准单点定位的精度通常只能达到米级。因此，标准单点定位主要用于导航定位及资源调查、勘测等低精度定位领域。利用精密星历和精密卫星钟差数据，由相位观测数据进行的定位方式，称为精密单点定位，精密单点定位的精度可以达到厘米量级。
在影响单点定位结果的三类误差源中，与信号传播有关的误差(主要是指电离层延迟和对流层延迟误差)无法很好地被模型化，因此，对定位的影响较大。
对于电离层延迟，单频接收机需要通过模型加以改正；双频接收机则可以通过线性组合法消除电离层延迟的影响；对于对流层延迟，由于对流层的性质为非弥散性，故不能用双频观测数据线性组合的方法消除对流层延迟误差，其对定位精度的影响不可以忽略。
对流层延迟分为干延迟和湿延迟两部分。在通常的定位解算策略中，天顶方向的干延迟由对流层延迟模型计算得出，通过映射函数将其映射到传播路径；而天顶方向的湿延迟则作为未知参数求解，并通过映射函数映射到传播路径。GNSS测量中对流层延迟模型或方法有萨斯塔莫宁模型、UNB3模型、气象元素法等。
以上三种对流层延迟改正方法存在如下问题：萨斯塔莫宁模型的精度相对较高，但是需要用到实测的气象数据，气象数据的获取问题，限制了萨斯塔莫宁模型的应用；UNB3模型不需要实测气象数据注入，但是UNB3模型的精度有限；气象元素法需要用到外部气象设备采集数据，不同设备的精度存在差异且设备笨重、昂贵，增加了外业采集数据的工作量；同时气象元素法也降低了单点定位的工作效率，很难为全天候的单点定位提供所需数据。
以上三种对流层延迟改正方法构建多种模型，以适应不同条件下改正对流层延迟的需要，考虑到实用性和普适性会使得测量外业和内业工作量增大，数据处理效率降低。此外，对于单点定位来说，对流层的延迟误差校正多采用经验模型，此类经验模型的参数估计也多采用 经验值，因而，难以符合实际单点定位中的参数估计。
发明内容
本发明的目的在于提出一种基于球谐展开的GNSS标准单点定位方法，以便能够快速、高效、准确地进行标准单点定位。
本发明为了实现上述目的，采用如下技术方案：
一种基于球谐展开的GNSS标准单点定位方法，采用球谐展开描述与测站和卫星之间的高度角和方位角有关的误差项，包括对流层延迟误差以及电离层延迟误差，使用GNSS观测数据和卫星星历数据即可实现测站空间位置的计算，进而实现实现GNSS标准单点定位。
本发明具有如下优点：
相比于现有的利用经验模型进行对流层延迟改正方法，本发明方法能够高效、快速地获得测点的位置信息，具有实际操作简单、数据处理便捷以及计算效率高等优势。
附图说明
图1为本发明实施例1中基于球谐展开的GNSS标准单点定位的流程示意图；
图2为本发明实施例2中基于球谐展开的GNSS精密单点定位的流程示意图；
图3为本发明实施例基于球谐展开的GNSS标准单点定位方法的结果统计图；
图4为本发明实施例基于球谐展开的GNSS精密单点定位方法的结果统计图。
具体实施方式
下面结合附图以及具体实施方式对本发明作进一步详细说明：
实施例1
本实施例述及了一种基于球谐展开的GNSS标准单点定位方法。如图1所示，该基于球谐展开的GNSS标准单点定位方法包括如下步骤：
I.1.建立利用伪距观测值进行的标准单点定位观测方程。
实现标准单点定位，读取观测数据和广播星历数据，处理伪距观测数据，计算与测站和卫星间的高度角和方位角无关或关系不密切的误差，包括卫星星历误差和卫星钟差。由于卫星星历误差的量级与标准单点定位误差的量级基本相同，因此，广播星历可用于低精度的标准单点定位，在标准单点定位中卫星星历误差可以暂时忽略。
在标准单点定位中，由广播星历内插计算得到每个历元对应的卫星位置；而卫星钟差使用广播星历提供的卫星钟差参数求得，一般卫星钟差t s由一个多项式拟合得到：t s＝α 0+α 1(t i-t 0)+α 2(t i-t 0) 2，α 0、α 1和α 2分别为钟差、钟漂和钟的频漂系数，可以在导航电文中读取；t 0为多项式计算卫星钟差的参考时刻，t i表示信号接收时刻即数据记录时刻。
地面接收机钟大多采用石英钟，地面接收机的精度低于卫星钟，钟差变化无规律性，无法通过多项式较好拟合，因此，将接收机钟差t r作为待求参数进行估计。
建立利用伪距观测值进行的标准单点定位观测方程，如公式(1)所示：
Figure PCTCN2021123845-appb-000001
式中，i表示第i观测历元，j表示卫星编号；ρ表示伪距观测值，
Figure PCTCN2021123845-appb-000002
表示卫星j在第i观测历元的伪距观测值；
Figure PCTCN2021123845-appb-000003
表示接收机的位置与卫星j在第i观测历元的位置之间的几何距离；c表示光速；t r表示测站接收机钟差，(t r) i表示第i观测历元测站的接收机钟差；t s表示卫星钟差，
Figure PCTCN2021123845-appb-000004
表示卫星j在第i观测历元的钟差；
Figure PCTCN2021123845-appb-000005
表示卫星j在第i观测历元的信号传播路径上的对流层延迟误差；
Figure PCTCN2021123845-appb-000006
表示卫星j在第i观测历元的信号传播路径上的电离层延迟误差；ε ρ表示伪距观测数据残差；定义卫星j在第i观测历元观测瞬间的空间三维坐标为
Figure PCTCN2021123845-appb-000007
测站的近似坐标为(X 0,Y 0,Z 0)，则卫星到测站近似位置的几何距离
Figure PCTCN2021123845-appb-000008
表示为：
Figure PCTCN2021123845-appb-000009
GNSS伪距观测数据类型分为单频观测数据和双频观测数据。若进行标准单点定位的GNSS接收机为双频接收机，则电离层延迟误差通过双频消电离层组合来消除，组合形式如下：
Figure PCTCN2021123845-appb-000010
ρ为伪距消电离层组合观测值；f 1、f 2为相应观测值的频率；ρ 1、ρ 2分别表示两种频率对应的伪距观测值。若观测数据为双频观测值时，则公式(1)中，电离层延迟误差的一阶主项通过双频消电离层组合消除，电离层延迟误差项
Figure PCTCN2021123845-appb-000011
为0。此时，本实施例1中的标准单点定位观测方程中，球谐展开主要用于表示对流层延迟误差项，球谐展开的阶次与单频伪距观测值构建的观测方程有所不同。
I.2.在标准单点定位中，对流层延迟误差和电离层延迟误差都是卫星信号穿过大气层时产生的一种信号延迟误差，都与卫星信号的传播有关；而卫星信号的传播路径又与测站和卫星的位置有关，利用测站坐标和卫星位置可以计算出测站和卫星之间的高度角和方位角，应用球谐展开，就能够表示与测站和卫星间的高度角和方位角有关的误差项，包括对流层延迟误差和电离层延迟误差，以消除或最大程度的削弱它们对标准单点定位结果的影响。基于球谐展开的标准单点定位观测方程表示为：
Figure PCTCN2021123845-appb-000012
式中，n为球谐展开的阶数，m为球谐展开的次数，Nmax为球谐展开的最大阶数；
Figure PCTCN2021123845-appb-000013
表示n阶m次的缔合勒让德多项式；C nm和S nm分别表示n阶m次的球谐函数模型的系数，C nm和S nm为球谐展开的待求参数；
Figure PCTCN2021123845-appb-000014
表示测站与卫星j在第i观测历元之间的高度角，
Figure PCTCN2021123845-appb-000015
表示测站与卫星j在第i观测历元之间的方位角。
经过公式(2)表示的基于球谐展开的标准单点定位观测方程，不涉及对流层延迟误差的改正模型，无需考虑不同模型在不同地区的适用问题，并且本实施例将球谐展开的系数当作待求参数直接求解，来改正与测站和卫星之间的高度角和方位角有关的误差项，主要包括对流层延迟误差，无需考虑气象信息，更加具有普适性。为了方便表示球谐函数展开，令：
Figure PCTCN2021123845-appb-000016
则公式(2)简化为：
Figure PCTCN2021123845-appb-000017
由简化后的标准单点定位观测方程(3)，得到标准单点定位误差方程(4)，公式如下：
Figure PCTCN2021123845-appb-000018
其中，v ρ表示伪距观测数据的改正值；
Figure PCTCN2021123845-appb-000019
表示卫星j在第i观测历元的伪距观测数据的改正值。
I.3.基于公式(4)得到标准单点定位误差方程的线性化表达式，并进行简化处理，对简化后的标准单点定位误差方程的线性化表达式进行滑动解算。
I.3.1.对标准单点定位误差方程(4)在测站的近似坐标(X 0,Y 0,Z 0)处进行泰勒级数展开并保留一阶项，则标准单点定位误差方程的线性化表达式为：
Figure PCTCN2021123845-appb-000020
公式(5)中，令
Figure PCTCN2021123845-appb-000021
Figure PCTCN2021123845-appb-000022
Figure PCTCN2021123845-appb-000023
Figure PCTCN2021123845-appb-000024
表示测站近似坐标与卫星j在第i观测历元坐标计算出的待求参数dX的系数，待求参数dX是测站近似坐标X 0的改正数；
Figure PCTCN2021123845-appb-000025
表示测站近似坐标与卫星j在第i观测历元坐标计算出的待求参数dY的系数，待求参数dY是测站近似坐标Y 0的改正数；
Figure PCTCN2021123845-appb-000026
表示测站近似坐标与卫星j在第i观测历元坐标计算出的待求参数dZ的系数，待求参数dZ是测站近似坐标Z 0的改正数；
Figure PCTCN2021123845-appb-000027
表示由测站近似坐标与卫星j在第i观测历元坐标计算出的球谐展开的待求参数C nm的系数，
Figure PCTCN2021123845-appb-000028
表示由测站近似坐标与卫星j在第i观测历元坐标计算出的球谐展开的待求参数S nm的系数；t r表示测站接收机钟差，因为测站接收机钟差无法很好地被多项式拟合，所以将其作为待求参数进行解算，作为待求参数求解时，dt r的系数为1；d(t r) i表示在第i观测历元测站接收机钟差的待求参数。
简化后的标准单点定位误差方程的线性化表达式，如公式(6)所示；
Figure PCTCN2021123845-appb-000029
公式(6)中，i＝1,2,…,epoch，epoch表示最大观测历元数；
Figure PCTCN2021123845-appb-000030
表示由卫星j在第i观测历元伪距观测值、卫星j在第i观测历元的钟差以及测站近似坐标与卫星j在第i观测历元之间的几何距离计算出的常数项，
Figure PCTCN2021123845-appb-000031
由公式(6)能够看出，此时，有epoch×j个观测方程。而待求参数包括3个位置参数(dX,dY,dZ)，epoch个接收机钟差dt r参数以及(Nmax+1) 2个球谐展开的参数C nm、S nm。
当epoch充分多时，能够通过最小二乘求得所有未知数的最优解。因此只要改正项是可以线性描述的，就能通过增加观测历元的方式来增加观测方程数量，从而解算出这些未知参数。
通过以上分析，为了确保未知参数的估计为最优解，本实施例选择合适的滑动窗口，对简化后的标准单点定位误差方程的线性表达式(6)进行滑动解算。
设定滑动窗口选择i个历元，每个历元最多能够观测到j颗卫星，则在这个滑动窗口下所有可视卫星的伪距观测值组成imax个方程，其中，imax＝i×j。
定义标准单点定位误差方程的系数矩阵、伪距观测值的改正数向量矩阵、常数项矩阵以及未知参数矩阵分别为矩阵B、V、L、X，则它们的表达式分别如下：
Figure PCTCN2021123845-appb-000032
Figure PCTCN2021123845-appb-000033
Figure PCTCN2021123845-appb-000034
X＝[dX dY dZ (dt r) 1 … (dt r) i C 00 … C NmaxNmax S NmaxNmax] T。
则标准单点定位误差方程的线性化表达式(6)用一个矩阵形式表达，如公式(7)所示：
V＝BX-L      (7)
在构建公式(7)的系数矩阵B时，只需要通过简单的数学计算，计算每一个待求参数(位置参数、测站接收机钟差参数及球谐展开的参数)的系数即可，计算过程简单，构建误差方程的速度迅速，实现标准单点定位的效率高。目前，与测站和卫星之间的高度角和方位角有关的误差项，如对流层延迟误差等还无法用模型精确改正，参考处理测站接收机钟差参数的方式，将其当做待求参数进行求解，是一种很好地方式。基于球谐展开的GNSS标准单点定位方法的提出正是基于这一思想，利用球谐展开表示与测站和卫星之间的高度角和方位角有关的误差项，在标准单点定位方程中直接求解球谐展开的系数。
未知参数X的最小二乘估值为：X＝(B TPB) -1B TPL    (8)
公式(8)中，P为单位权矩阵。
I.3.2.如果标准单点定位误差方程的系数矩阵B是良态，则用公式(8)求出的未知参数的解是最优无偏估计，当系数矩阵B是病态时，则利用最小二乘法求出的不再是最优解。
因此，需要首先判断标准单点定位误差方程的系数矩阵B是否为病态，若系数矩阵B是病态，则执行步骤I.3.3；否则，若系数矩阵B为良态，执行步骤I.3.4。
本实施例通过计算的系数阵B的条件数，设置条件数的经验阈值，并将条件数与经验阈值作比较，当条件数大于经验阈值时，则表明误差方程(7)为病态，否则为良态。
I.3.3.为了解决误差方程的病态问题，本实施例提出了利用截断奇异值分解的方法为最小 二乘谱修正迭代提供待求参数X的初值，基于最小二乘谱修正迭代计算待求参数X的值。截断奇异值法的思想是计算矩阵S中奇异值的均值，将该均值作为截断奇异值的阈值；筛选出小于阈值的奇异值将其零化。由于截断奇异值分解方法是将方程的最小二乘法转换为直接解法，因而避免了病态方程的解过分偏离真实解，有效解决了误差方程系数矩阵B的病态问题。
设系数矩阵B∈R n×m，R n×m表示n行m列的实数阵；则系数矩阵B的奇异值分解为：
B＝USV T    (9)
公式(9)中，U∈n×n，V∈m×m，U、V均为正交矩阵，S∈n×m为一对角矩阵。
系数矩阵B∈R n×m的截断奇异值矩阵B k定义为：
Figure PCTCN2021123845-appb-000035
公式(10)中，矩阵S中最小的(r-k)个非奇异值用零代替，即被截断，且k≤r；r表示系数矩阵B的秩，k表示矩阵S中保留的奇异值数量；u i表示矩阵U对应的向量，v i表示矩阵V对应的向量；σ i表示矩阵S中保留的奇异值。计算矩阵S中奇异值的均值，将该值作为截断奇异值的阈值；矩阵S中奇异值大于该阈值的保留，小于该阈值的进行零化处理，则公式(7)的截断奇异值解
Figure PCTCN2021123845-appb-000036
为：
Figure PCTCN2021123845-appb-000037
公式(11)中，
Figure PCTCN2021123845-appb-000038
根据最小二乘原理V TPV＝min，公式(8)写成以下形式：
(B TPB)X＝(B TPL)    (12)
将公式(12)的等式左端和右端同时加未知参数KX，并化简得到公式(13)：
(B TPB+KI)X＝B TPL+KX    (13)
由公式(13)整理得，最小二乘谱修正迭代公式为：
X＝(B TPB+KI) -1(B TPL+KX)    (14)
公式(14)中，K为任意实数。最小二乘谱修正迭代是一种基于最小二乘的迭代计算方法；基于最小二乘谱修正迭代计算未知参数X的解。求解过程如下：
基于最小二乘谱修正迭代解算X包括迭代①和迭代②两个迭代过程；其中，设定用于判断迭代①是否收敛的第一比较阈值，设定用于判断迭代②是否收敛的第二比较阈值。
迭代①第一次迭代时，设定K的初始值为1，将公式(11)得到的截断奇异值解
Figure PCTCN2021123845-appb-000039
作为公式(14)第一次迭代时X的初值，并赋值给公式(14)等式右边的X；
每次迭代过程中，将上一次迭代时利用公式(14)得到的未知参数X的值，作为本次迭代过程中X的初值，并赋值给公式(14)等式右边的X；
基于公式(14)求解得到每次迭代过程中未知参数X的值；
将每次迭代过程中求解得到的X的值，与每次迭代过程中X的初值进行差值运算；
若经过差值运算，得到的差值大于第一比较阈值，则表示迭代不收敛，需要改正最小二乘谱修正迭代中的K值，即将K值递增1，继续上述迭代①过程；
若经过差值运算得到的差值小于或等于第一比较阈值，则跳出迭代①，进入迭代②过程。
迭代②将利用公式(14)求解得到的未知参数X的值与第二比较阈值进行比较；
若利用公式(14)求解得到的未知参数X的值大于第二比较阈值，则表示迭代②未收敛，此时，将位置参数(dX,dY,dZ)分别加到测站的近似坐标(X 0,Y 0,Z 0)上，更新测站的近似坐标，通过公式(5)对更新测站近似坐标后的标准单点定位误差方程进行线性化，通过公式(6)简化线性化后的标准单点定位误差方程；重新进入迭代①过程求解出未知参数X的值。
若利用公式(14)求解得到的未知参数X的值小于或等于第二比较阈值，则表示迭代②收敛，跳出迭代；此时利用公式(14)得到的未知参数X的值，即未知参数X的解。
求得未知参数X的解后，转到步骤I.3.5。
I.3.4.利用最小二乘法求解未知参数X的解，转到步骤I.3.5。
步骤I.3.4中，基于最小二乘法求解未知参数X的解的过程如下：
设定用于判断迭代是否收敛的第三比较阈值；每次迭代过程中，利用公式(8)求解出未知参数X的值，将本次迭代得到的未知参数X的值与第三比较阈值进行比较；
若利用公式(8)求解得到的未知参数X的值大于第三比较阈值，则表示迭代未收敛，此时，将位置参数(dX,dY,dZ)，分别加到测站的近似坐标(X 0,Y 0,Z 0)上，更新测站的近似坐标，执行公式(5)-公式(8)的计算过程，并重新求解未知参数X的值；
若利用公式(8)求解得到的未知参数X的值小于或等于第三比较阈值，则表示迭代收敛，跳出迭代；此时，利用公式(8)求解得到的未知参数X的值，即未知参数X的解。
I.3.5.根据计算得到的未知参数X的解，得到未知参数X中包含的参数值，即测站的位置参数(dX,dY,dZ)、接收机钟差dt r以及球谐展开的系数C nm和S nm；
其中，未知参数X中的位置参数(dX,dY,dZ)分别加到测站的近似坐标(X 0,Y 0,Z 0)上，即 标准单点定位求出的地心地固坐标系下的测站坐标(X,Y,Z)。
通过求解得到的球谐展开的系数C nm和S nm可直接用于后续标准单点定位。
具体应用过程为：读取已有的球谐展开系数，通过球谐展开表示与测站和卫星之间的高度角和方位角有关的误差项，主要包括对流层延迟误差和电离层延迟误差；构建单点定位观测方程时，不再考虑与测站和卫星之间的高度角和方位角有关的误差项的改正问题，省略定位中对流层延迟误差和电离层延迟误差的求解过程，使得定位方程仅仅需要解算测站的位置参数、测站接收机钟差参数即可，极大地简化了标准单点定位流程。
I.3.6.将测站坐标由地心地固坐标系转换到站心坐标系，实现GNSS标准单点定位。
将测站坐标由地心地固坐标系转换到站心坐标系的过程如下：
地心地固坐标系转经纬高坐标系的公式为：
Figure PCTCN2021123845-appb-000040
式中，(X,Y,Z)为在地心地固坐标系的坐标，(lon,lat,alt)为在经纬高坐标系的坐标；M表示基准椭球体的曲率半径，E表示椭球偏心率；根据测站近似坐标(X 0,Y 0,Z 0)，利用公式(15)计算得到该测站近似坐标(X 0,Y 0,Z 0)在经纬高坐标系下的对应坐标，即(lon 0,lat 0,alt 0)。
地心地固坐标系转换到站心坐标系的公式为：
Figure PCTCN2021123845-appb-000041
式中，(e 1,n 1,u 1)为以测站近似坐标(X 0,Y 0,Z 0)为原点求出的(X,Y,Z)坐标的位置。
S为坐标转换矩阵，
Figure PCTCN2021123845-appb-000042
根据以上坐标转换关系，将测站坐标由地心地固坐标系转换到站心坐标系。
本实施例1中的GNSS标准单点定位方法，基于球谐展开表示与测站和卫星间的方位角和高度角有关的误差项，从而能够快速、高效、准确地进行标准单点定位。
实施例2
本实施例2述及了一种基于球谐展开的GNSS精密单点定位方法。如图2所示，该基于球谐展开的GNSS精密单点定位方法包括如下步骤：
II.1.建立利用载波相位观测值进行的精密单点定位观测方程。
读取观测数据、广播星历数据以及精密星历数据，对载波相位观测数据进行预处理，计算与测站和卫星间的高度角和方位角无关或关系不密切的误差；以上误差主要包括计算卫星钟差、探测周跳以及计算载波相位整周模糊度参数的初始值等。
对于精密单点定位，广播星历的精度无法满足精密单点定位的要求，因此必须采用卫星精密星历计算卫星位置；利用二阶多项式估计的卫星钟差无法满足精密单点定位的精度要求；由于不同的卫星钟差均不相同，必须事先确定其大小，再代入到观测方程消除其影响。
目前，IGS提供的15min、5min采样间隔的事后精密星历和5min、30s采样间隔的卫星钟差产品，能够满足精密定位的精度要求。观测值的采样间隔一般都小于上述文件中数据的采样间隔，因此，需要内插计算得到每个历元对应的卫星位置和卫星钟差。
地面接收机钟大多采用石英钟，精度低于卫星钟，钟差变化无规律性，无法通过多项式较好拟合，因此，将接收机钟差t r作为待求参数进行估计。
建立的利用载波相位观测值进行的精密单点定位观测方程，如公式(1)所示。
Figure PCTCN2021123845-appb-000043
公式(1)中，i表示第i观测历元；j表示卫星编号；
Figure PCTCN2021123845-appb-000044
表示载波相位观测值，将载波相位原始观测值乘上相应的波长λ，得到以距离表示的载波相位观测数据，称为等效距离；
Figure PCTCN2021123845-appb-000045
表示卫星j在第i观测历元的载波相位观测值，即等效距离；
Figure PCTCN2021123845-appb-000046
表示接收机的位置与卫星j在第i观测历元位置之间的几何距离；t r表示接收机钟差，(t r) i表示第i观测历元接收机钟差；t s表示卫星钟差，
Figure PCTCN2021123845-appb-000047
表示卫星j在第i观测历元的钟差；
Figure PCTCN2021123845-appb-000048
表示卫星j在第i观测历元的信号传播路径上的对流层延迟误差；
Figure PCTCN2021123845-appb-000049
表示卫星j在第i观测历元的信号传播路径上的电离层延迟误差；c表示光速；λ表示载波的波长；N表示卫星相位观测数据的整周模糊度参数，
Figure PCTCN2021123845-appb-000050
表示卫星j在第i观测历元的相位观测数据整周模糊度参数；
Figure PCTCN2021123845-appb-000051
表示载波相位观测数据残差。定义卫星j在第i历元观测瞬间的空间三维坐标为
Figure PCTCN2021123845-appb-000052
测站的近似坐标为(X 0,Y 0,Z 0)，则卫星到测站近似位置的几何距离
Figure PCTCN2021123845-appb-000053
表示为：
Figure PCTCN2021123845-appb-000054
GNSS载波相位观测数据类型分为单频观测数据和双频观测数据。
若进行标准单点定位的GNSS接收机为双频接收机，则电离层延迟误差可以通过双频消电离层组合来消除，具体的组合形式如下：
Figure PCTCN2021123845-appb-000055
式中，
Figure PCTCN2021123845-appb-000056
表示载波相位消电离层组合观测值(等效距离)；f 1、f 2表示相应观测值的频率；
Figure PCTCN2021123845-appb-000057
分别表示两种频率对应的载波相位观测值(等效距离)。若观测数据为双频观测值时，则上述公式(1)中，电离层延迟误差的一阶主项通过双频消电离层组合消除，电离层延迟误差项
Figure PCTCN2021123845-appb-000058
为0。
此时，本实施例2中的精密单点定位观测方程中，球谐展开主要用于表示对流层延迟误差项，球谐展开的阶次与单频载波相位观测值构建的观测方程有所不同。
在精密单点定位中，利用载波相位观测值组成GF(几何无关)组合观测值，进行整周模糊度的解算和周跳的探测。具体的组合形式如下：
Figure PCTCN2021123845-appb-000059
式中，
Figure PCTCN2021123845-appb-000060
表示GF组合载波相位观测值(等效距离)；
Figure PCTCN2021123845-appb-000061
分别表示两种频率对应的载波相位观测值；λ 1、λ 2分别表示两种频率对应的载波相位的波长；f 1、f 2表示两种载波相位的频率。GF组合观测值与几何距离无关，因此，该组合适合用来探测周跳。
II.2.基于球谐展开表示与测站和卫星间的高度角和方位角有关的误差项，其中，误差项主要包括对流层延迟误差以及电离层延迟误差。在精密单点定位中，对流层延迟误差和电离层延迟误差都是卫星信号穿过大气层时产生的一种信号延迟误差，都与卫星信号的传播有关；而卫星信号的传播路径又与测站和卫星的位置有关，利用测站坐标和卫星位置可以计算出测站和卫星之间的高度角和方位角，应用球谐展开，就可以表示与测站和卫星间的高度角和方位角有关的误差项，主要包括对流层延迟误差和电离层延迟误差，以消除或最大程度的削弱其对精密单点定位结果的影响。经过球谐展开后的精密单点定位观测方程表示为：
Figure PCTCN2021123845-appb-000062
式(2)中，n为球谐展开的阶数，m为球谐展开的次数，Nmax为球谐展开的最大阶数；P nm(cose)表示n阶m次的缔合勒让德多项式；C nm和S nm表示n阶m次的球谐展开的系数；
Figure PCTCN2021123845-appb-000063
表示测站与卫星j在第i观测历元之间的高度角，
Figure PCTCN2021123845-appb-000064
表示测站与卫星j在第i观测历元之间的方位角。公式(2)提出的基于球谐展开的精密单点定位观测方程，不涉及对流层延迟误差的 改正模型，无需考虑不同模型在不同地区的适用问题，并且本实施例将球谐展开的系数当作待求参数直接求解，来改正与测站和卫星之间的高度角和方位角有关的误差项，主要包括对流层延迟误差，无需考虑气象信息，更加具有普适性；为了方便表示球谐函数展开，令：
Figure PCTCN2021123845-appb-000065
则公式(2)简化为：
Figure PCTCN2021123845-appb-000066
由公式(3)中的精密单点定位观测方程，得到精密单点定位误差方程：
Figure PCTCN2021123845-appb-000067
其中，
Figure PCTCN2021123845-appb-000068
表示载波相位观测数据(等效距离)的改正值。
Figure PCTCN2021123845-appb-000069
表示第i观测历元卫星j的载波相位观测数据(等效距离)的改正值。
II.3.基于公式(4)得到精密单点定位误差方程的线性化表达式，并简化处理，然后对简化后的精密单点定位误差方程的线性化表达式进行滑动解算。
II.3.1.对精密单点定位误差方程(4)在测站的近似坐标(X 0,Y 0,Z 0)处进行泰勒级数展开并保留一阶项，则精密单点定位误差方程的线性化表达式为：
Figure PCTCN2021123845-appb-000070
公式(5)中，
Figure PCTCN2021123845-appb-000071
表示卫星j在第i历元的模糊度初值；
Figure PCTCN2021123845-appb-000072
表示卫星j在第i历元的相位观测数据整周模糊度参数，dN表示相位观测数据整周模糊度参数的改正数，
Figure PCTCN2021123845-appb-000073
表示卫星j在第i历元的相位观测数据整周模糊度参数的改正数，在精密单点定位中，由于整周模糊度参数不能被很好地模型化，因此将其改正数作为待求参数进行求解，作为待求参数求解时，
Figure PCTCN2021123845-appb-000074
的系数为1；d(t r) i表示在第i观测历元测站接收机钟差的待求参数，与整周模糊度的处理方式相似，因为测站接收机钟差无法很好地被多项式拟合，因此将其作为待求参数进行求解，作为待求参数求解时，dt r的系数为1。
令
Figure PCTCN2021123845-appb-000075
Figure PCTCN2021123845-appb-000076
Figure PCTCN2021123845-appb-000077
Figure PCTCN2021123845-appb-000078
为测站近似坐标与卫星j在第i观测历元坐标计算出的待求参数dX的系数，待求参数dX是测站近似坐标X 0的改正数；
Figure PCTCN2021123845-appb-000079
为测站近似坐标与卫星j在第i观测历元坐标计算出的待求参数dY的系数，待求参数dY是测站近似坐标Y 0的改正数；
Figure PCTCN2021123845-appb-000080
为测站近似坐标与卫星j在第i观测历元坐标计算出的待求参数dZ的系数，待求参数dZ是测站近似坐标Z 0的改正数；A nm表示球谐展开待求参数C nm的系数，B nm表示球谐展开待求参数S nm的系数；
Figure PCTCN2021123845-appb-000081
表示由测站近似坐标与卫星j在第i观测历元坐标计算出的球谐展开的待求参数C nm的系数，
Figure PCTCN2021123845-appb-000082
为由测站近似坐标与卫星j在第i观测历元坐标计算出的球谐展开的待求参数S nm的系数。
简化后的精密单点定位误差方程的线性化表达式，如公式(6)所示：
Figure PCTCN2021123845-appb-000083
公式(6)中，i＝1,2,,…,epoch，epoch表示最大观测历元；
Figure PCTCN2021123845-appb-000084
表示由卫星j在第i观测历元载波相位观测值、卫星j在第i观测历元的钟差、卫星j在第i历元的相位观测数据整周模糊度初始值以及测站近似坐标与卫星j在第i观测历元之间的几何距离计算出的常数项，
Figure PCTCN2021123845-appb-000085
由公式(6)可以看出，此时有epoch×j个观测方程。而待求参数包括：3个位置参数(dX,dY,dZ)，epoch个接收机钟差参数dt r，j个模糊度参数dN以及(Nmax+1) 2个球谐展开的系数C nm、S nm。
当epoch充分多时，能够通过最小二乘求得所有未知数的最优解，只要改正项是可以线性描述的，就可以通过增加观测历元的方式来增加观测方程数量，从而解算出这些未知参数。
设定滑动窗口选择i个历元，每个历元最多能够观测到j颗卫星，则在这个滑动窗口下所有可视卫星的载波相位观测值组成imax个方程，其中，imax＝i×j。
定义精密单点定位误差方程的系数矩阵、载波相位观测值的改正数向量矩阵、常数项矩阵以及未知参数矩阵分别为矩阵B、V、L、X，则它们的表达式分别如下：
Figure PCTCN2021123845-appb-000086
Figure PCTCN2021123845-appb-000087
Figure PCTCN2021123845-appb-000088
Figure PCTCN2021123845-appb-000089
Figure PCTCN2021123845-appb-000090
表示第1颗卫星在第1个观测历元的整周模糊度参数，
Figure PCTCN2021123845-appb-000091
表示第j颗卫星在第i个观测历元的整周模糊度参数，精密单点定位中求得的是整周模糊度初始值的改正数。
则精密单点定位误差方程的线性化表达式(6)用一个矩阵形式表达，如公式(7)所示；
V＝BX-L   (7)
构建精密单点定位误差方程式(7)的系数矩阵B时，只需要通过简单的数学计算，计算每一个待求参数(位置参数、测站接收机钟差参数、整周模糊度参数及球谐展开的参数)的系数即可，计算过程简单，构建误差方程的速度迅速，实现精密单点定位的效率高。目前，与测站和卫星之间的高度角和方位角有关的误差项，如对流层延迟误差等还无法用模型精确改正，参考处理测站接收机钟差参数的方式，将其当做待求参数进行求解，是一种很好地方式；基于球谐展开的GNSS精密单点定位方法的提出正是基于这一思想，利用球谐展开表示与测站和卫星之间的高度角和方位角有关的误差项，在精密单点定位方程中直接求解球谐展开的系数。未知参数X的最小二乘估值为：X＝(B TPB) -1B TPL   (8)
公式(8)中，P为单位权矩阵。
II.3.2.如果精密单点定位误差方程的系数矩阵B是良态，则用公式(8)求出的未知参数的解时最优无偏估计，当系数矩阵B是病态时，则利用最小二乘法求出的不再是最优解。
因此，需要首先判断精密单点定位误差方程的系数矩阵B是否为病态，若系数矩阵B是病态，则执行步骤II.3.3；否则，系数矩阵B为良态，执行步骤II.3.4。
本实施例通过计算精密单点定位误差方程的系数阵B的条件数，设置条件数的经验阈值，将条件数与经验阈值作比较，当条件数大于经验阈值时，则误差方程为病态。
II.3.3.为了解决误差方程的病态问题，本实施例提出了利用截断奇异值分解的方法，为最小二乘谱修正迭代提供待求参数X的初值，基于最小二乘谱修正迭代计算待求参数X的 值。由于截断奇异值分解方法是将方程的最小二乘法转换为直接解法，因而避免了病态方程的解过分偏离真实解，有效解决了误差方程系数矩阵B的病态问题。
设系数矩阵B∈R n×m，R n×m表示n行m列的实数阵，则系数矩阵B的奇异值分解为：
B＝USV T   (9)
公式(9)中，U∈n×n，V∈m×m，U、V均为正交矩阵，S∈n×m为一对角矩阵。
系数矩阵B∈R n×m的截断奇异值矩阵B k定义为：
Figure PCTCN2021123845-appb-000092
公式(10)中，矩阵S中最小的(r-k)个非奇异值用零代替，即被截断，且k≤r；r表示系数矩阵B的秩，k表示矩阵U中保留的奇异值数量，u i表示矩阵U对应的向量，v i表示矩阵V对应的向量；σ i表示矩阵U中保留的奇异值。计算矩阵S中奇异值的均值，将该值作为截断奇异值的阈值；矩阵S中奇异值大于该截断奇异值的阈值的保留，小于该截断奇异值的阈值的进行零化处理。则公式(7)的截断奇异值解
Figure PCTCN2021123845-appb-000093
为：
Figure PCTCN2021123845-appb-000094
公式(11)中，
Figure PCTCN2021123845-appb-000095
最小二乘谱修正迭代是一种基于最小二乘的迭代计算方法。基于最小二乘谱修正迭代计算待求参数X；根据最小二乘原理V TPV＝min，公式(8)写成以下形式：
(B TPB)X＝(B TPL)    (12)
将公式(12)的等式左端和右端同时加未知参数KX，并化简得到公式(13)：
(B TPB+KI)X＝B TPL+KX    (13)
由公式(13)整理得，最小二乘谱修正迭代公式为：
X＝(B TPB+KI) -1(B TPL+KX)    (14)
公式(14)中，K为任意实数；基于最小二乘谱修正迭代计算未知参数X的解。
基于最小二乘谱修正迭代公式计算X的解的过程如下：
基于最小二乘谱修正迭代解算X包括迭代①和迭代②两个迭代过程；其中，设定用于判断迭代①是否收敛的第一比较阈值，设定用于判断迭代②是否收敛的第二比较阈值；
迭代①第一次迭代时，设定K的初始值为1，将公式(11)得到的截断奇异值解
Figure PCTCN2021123845-appb-000096
作为公式(14)第一次迭代时X的初值，并赋值给公式(14)等式右边的X；
每次迭代过程中，将上一次迭代时利用公式(14)得到的未知参数X的值，作为本次迭代过程中X的初值，并赋值给公式(14)等式右边的X；
基于公式(14)求解得到每次迭代过程中未知参数X的值；
将每次迭代过程中求解得到的X的值，与每次迭代过程中X的初值进行差值运算；
若经过差值运算，得到的差值大于第一比较阈值，则表示迭代不收敛，需要改正最小二乘谱修正迭代中的K值，即将K值递增1，继续上述迭代①过程；
若经过差值运算得到的差值小于或等于第一比较阈值，则跳出迭代①，进入迭代②过程；
迭代②将利用公式(14)求解得到的未知参数X的值与第二比较阈值进行比较；
若利用公式(14)求解得到的未知参数X的值大于第二比较阈值，则表示迭代②未收敛，此时，将位置参数(dX,dY,dZ)分别加到测站的近似坐标(X 0,Y 0,Z 0)上，更新测站近似坐标，通过公式(5)对更新测站近似坐标后的精密单点定位误差方程进行线性化，通过公式(6)简化线性化后的精密单点定位误差方程；重新进入迭代①过程求解出未知参数X的值；
若利用公式(14)求解得到的未知参数X的值小于或等于第二比较阈值，则表示迭代②收敛，跳出迭代；此时利用公式(14)得到的未知参数X的值，即待求解的未知参数X的解。
未知参数X的解后，转到步骤II.3.5。
II.3.4.利用最小二乘法求解未知参数X的解，转到步骤II.3.5。
步骤II.3.4中，基于最小二乘法求解未知参数X的解的过程如下：
设定用于判断迭代是否收敛的第三比较阈值；每次迭代过程中，利用公式(8)求解出未知参数X的值，将本次迭代得到的未知参数X的值与第三比较阈值进行比较；
若利用公式(8)求解得到的未知参数X的值大于第三比较阈值，则表示迭代未收敛，此时，将位置参数(dX,dY,dZ)分别加到测站的近似坐标(X 0,Y 0,Z 0)上，更新测站近似坐标；
执行公式(5)-公式(8)的计算过程，并重新求解未知参数X的值；
若利用公式(8)求解得到的未知参数X的值小于或等于第三比较阈值，则表示迭代收敛，跳出迭代；此时，利用公式(8)求解得到的未知参数X的值，即待求解的未知参数X的解。
II.3.5.根据计算得到的未知参数X的解，得到未知参数X中包含的参数值，即测站的位置参数(dX,dY,dZ)、接收机钟差dt r、模糊度参数dN以及球谐展开的系数C nm和S nm；
其中，位置参数(dX,dY,dZ)分别加到测站的近似坐标(X 0,Y 0,Z 0)上，即精密单点定位求 出的地心地固坐标系下的测站坐标(X,Y,Z)。通过求解得到的球谐展开的系数C nm和S nm，可以直接用于后续精密单点定位，具体应用过程如下：
读取已有的球谐展开系数，通过球谐展开表示与测站和卫星之间的高度角和方位角有关的误差项，主要包括对流层延迟误差和电离层延迟误差，构建的精密单点定位观测方程可以不再考虑与测站和卫星之间的高度角和方位角有关的误差项的改正问题，省略定位中的对流层延迟误差和电离层延迟误差的求解过程，这将使得定位方程仅仅需要解算测站的位置参数、测站接收机钟差参数以及整周模糊度参数即可，极大地简化了精密单点定位流程。
II.3.6.将测站坐标由地心地固坐标系转换到站心坐标系，实现GNSS精密单点定位。其中，本实施例2中坐标转换过程与上述实施例1中坐标转换过程相同，此处不再赘述。
本实施例2中述及的GNSS精密单点定位方法，基于球谐展开表示与测站和卫星间的方位角和高度角有关的误差项，从而能够快速、高效、准确地进行精密单点定位。
下面对利用本发明中基于球谐展开的GNSS单点定位方法所能达到的精度进行说明：
以上海佘山站为例：
获取2019年1月7日佘山站的GNSS观测数据和卫星星历数据。
GNSS观测数据采样间隔为30s，观测数据选择GPS卫星的双频伪距观测值以及双频载波相位观测值，站点的参考坐标为IGS提供的站点坐标。
1.首先利用本发明中基于球谐展开的GNSS标准单点定位方法获取站点位置信息。
观测数据选择GPS卫星双频伪距观测值，则电离层延迟误差通过伪距线性组合消除。
通过球谐展开表示与卫星和测站间的高度角和方位角有关的误差项，主要是对流层延迟误差。球谐展开各阶次的系数、以及测站的接收机钟差需要当作待求参数求解。
将球谐函数展开到3阶，滑动窗口选择10个历元。待求参数包括：点位坐标3个方向的改正值，10个接收机钟差(每个历元一个接收机钟差参数)，球谐展开到3阶3次的16个系数：A 00、A 10、A 11、B 11、A 20、A 21、B 21、A 22、B 22、A 30、A 31、B 31、A 32、B 32、A 33、B 33。
卫星的截止高度角设置为5°，最小二乘谱修正迭代的系数k设为2。
GNSS标准单点定位方法的定位结果如图3所示，统计结果见表1。
表1 标准单点定位精度(单位/m)
  E N U
MAX 1.2662 1.7738 4.9711
MIN -1.3027 -2.0436 -5.0853
MEAN 0.1101 -0.0136 0.3269
STD 0.3727 0.3952 1.3887
RMS 0.3886 0.3954 1.4267
表2 Bernese定位精度(单位/m)
  E N U
MAX 6.0917 5.7515 9.7484
MIN -3.8416 -5.8132 -11.1756
MEAN 0.1268 -0.1259 0.5520
STD 1.2744 1.6002 3.0886
RMS 1.2807 1.6052 3.1375
通过表1统计结果不难发现，本发明标准单点定位方法得到的定位结果在E方向上的RMS优于0.39m，在N方向上的RMS优于0.40m，在U方向上的RMS优于1.43m。
由瑞士伯尔尼大学天文研究所开发的国际知名GNSS高精度数据处理软件Bernese软件解算的标准单点定位结果，如表2所示。
可见，GNSS高精度数据处理软件Bernese软件解算的定位结果在E方向上的RMS优于1.28m，在N方向上的RMS优于1.61m，在U方向上的RMS优于3.14m。
从统计表中看出，本发明标准单点定位结果与Bernese软件标准单点定位结果精度相当。
2.利用本发明中基于球谐展开的GNSS标准单点定位方法获取站点位置信息。
观测数据选择GPS卫星双频载波相位观测值，则电离层延迟误差通过载波相位线性组合消除。通过球谐展开表示与卫星和测站间的高度角和方位角有关的误差项，主要是对流层延迟误差；球谐展开各阶次的系数需要当作待求参数求解。
每颗卫星载波相位观测数据数据的整周模糊度以及测站的接收机钟差也需要作为待求参数求解。将球谐函数展开到6阶，滑动窗口选择15个历元。
待求参数包括：点位坐标3个方向的改正值，15个接收机钟差(每个历元一个接收机钟差参数)，每个历元下每颗卫星的模糊度参数，球谐展开到6阶6次的49个系数：A 00、A 10、A 11、B 11、A 20、A 21、B 21、A 22、B 22、A 30、A 31、B 31、A 32、B 32、A 33、B 33、A 40、A 41、B 41、A 42、B 42、A 43、B 43、A 44、B 44、A 50、A 51、B 51、A 52、B 52、A 53、B 53、A 54、B 54、A 55、B 55、A 60、A 61、B 61、A 62、B 62、A 63、B 63、A 64、B 64、A 65、B 65、A 66、B 66。
卫星的截止高度角设置为5°，最小二乘谱修正迭代的系数k设为3。
GNSS精密单点定位方法的结果如图4所示，统计结果见表3。
表3 精密单点定位精度(单位/m)
  E N U
MAX 0.0230 0.0345 0.0432
MIN -0.0041 0.0173 0.0218
MEAN 0.0120 0.0246 0.0344
STD 0.0044 0.0031 0.0051
RMS 0.0128 0.0248 0.0348
表4 Bernese定位精度(单位/m)
  E N U
MAX 0.0274 0.0301 0.0461
MIN -0.0178 -0.0209 -0.0738
MEAN 0.0036 0.0048 -0.0082
STD 0.0075 0.0092 0.0200
RMS 0.0083 0.0104 0.0216
通过表3统计结果不难发现，本发明方法得到的定位结果在E方向上的RMS优于0.02m，在N方向上的RMS优于0.03m，在U方向上的RMS优于0.04m。由瑞士伯尔尼大学天文研究所开发的国际知名GNSS高精度数据处理软件Bernese软件解算的标准单点定位结果如表4所示。可见，GNSS高精度数据处理软件Bernese软件解算的定位结果在E方向上的RMS优于0.009m，在N方向上的RMS优于0.02m，在U方向上的RMS优于0.03m。
从统计表中看出，本发明精密单点定位结果与Bernese软件精密单点定位结果精度相当。
本发明基于球谐展开的GNSS单点定位方法，利用观测数据和卫星星历数据高效快速获得测点的点位信息，相比于现有的利用经验模型进行对流层延迟改正方法，具有解算效率高，无需提供气象文件，操作简便、工作难度小，测量结果可靠，测量精度高等优势。
当然，以上说明仅仅为本发明的较佳实施例，本发明并不限于列举上述实施例，应当说明的是，任何熟悉本领域的技术人员在本说明书的教导下，所做出的所有等同替代、明显变形形式，均落在本说明书的实质范围之内，理应受到本发明的保护。

Claims (4)

  1. 一种基于球谐展开的GNSS标准单点定位方法，其特征在于，包括如下步骤：
    I.1.建立利用伪距观测值进行的标准单点定位观测方程，如公式(1)所示：
    Figure PCTCN2021123845-appb-100001
    式中，i表示第i观测历元，j表示卫星编号；ρ表示伪距观测值，
    Figure PCTCN2021123845-appb-100002
    表示卫星j在第i观测历元的伪距观测值；
    Figure PCTCN2021123845-appb-100003
    表示接收机的位置与卫星j在第i观测历元的位置之间的几何距离；c表示光速；t r表示测站接收机钟差，(t r) i表示第i观测历元测站的接收机钟差；t s表示卫星钟差，
    Figure PCTCN2021123845-appb-100004
    表示卫星j在第i观测历元的钟差；
    Figure PCTCN2021123845-appb-100005
    表示卫星j在第i观测历元的信号传播路径上的对流层延迟误差；
    Figure PCTCN2021123845-appb-100006
    表示卫星j在第i观测历元的信号传播路径上的电离层延迟误差；ε ρ表示伪距观测数据残差；
    定义卫星j在第i观测历元观测瞬间的空间三维坐标为
    Figure PCTCN2021123845-appb-100007
    测站的近似坐标为(X 0,Y 0,Z 0)，则卫星到测站近似位置的几何距离
    Figure PCTCN2021123845-appb-100008
    表示为：
    Figure PCTCN2021123845-appb-100009
    I.2.基于球谐展开表示与测站和卫星之间的高度角和方位角有关的误差项，误差项包括对流层延迟误差以及电离层延迟误差；基于球谐展开的标准单点定位观测方程表示为：
    Figure PCTCN2021123845-appb-100010
    式中，n为球谐展开的阶数，m为球谐展开的次数，Nmax为球谐展开的最大阶数；
    Figure PCTCN2021123845-appb-100011
    表示n阶m次的缔合勒让德多项式；C nm和S nm分别表示n阶m次的球谐展开的系数，C nm和S nm为球谐展开的待求参数；
    Figure PCTCN2021123845-appb-100012
    分别表示测站与卫星j在第i观测历元之间的高度角和方位角；为了方便表示球谐展开，令：
    Figure PCTCN2021123845-appb-100013
    则公式(2)简化为：
    Figure PCTCN2021123845-appb-100014
    由简化后的标准单点定位观测方程(3)，得到标准单点定位误差方程(4)，公式如下：
    Figure PCTCN2021123845-appb-100015
    其中，v ρ表示伪距观测数据的改正值；
    Figure PCTCN2021123845-appb-100016
    表示卫星j在第i观测历元的伪距观测数据的改正值；
    I.3.基于公式(4)得到标准单点定位误差方程的线性化表达式，并进行简化处理，然后对简化后的标准单点定位误差方程的线性化表达式进行滑动解算；
    I.3.1.对标准单点定位误差方程(4)在测站的近似坐标(X 0,Y 0,Z 0)处进行泰勒级数展开并保留一阶项，则标准单点定位误差方程的线性化表达式如公式(5)所示；
    Figure PCTCN2021123845-appb-100017
    公式(5)中，令：
    Figure PCTCN2021123845-appb-100018
    Figure PCTCN2021123845-appb-100019
    Figure PCTCN2021123845-appb-100020
    其中，
    Figure PCTCN2021123845-appb-100021
    表示测站近似坐标与卫星j在第i观测历元坐标计算出的待求参数dX的系数，待求参数dX是测站近似坐标X 0的改正数；
    Figure PCTCN2021123845-appb-100022
    表示测站近似坐标与卫星j在第i观测历元坐标计算出的待求参数dY的系数，待求参数dY是测站近似坐标Y 0的改正数；
    Figure PCTCN2021123845-appb-100023
    表示测站近似坐标与卫星j在第i观测历元坐标计算出的待求参数dZ的系数，待求参数dZ是测站近似坐标Z 0的改正数；
    Figure PCTCN2021123845-appb-100024
    表示由测站近似坐标与卫星j在第i观测历元坐标计算出的球谐展开的待求参数C nm的系数；
    Figure PCTCN2021123845-appb-100025
    表示由测站近似坐标与卫星j在第i观测历元坐标计算出的球谐展开的待求参数S nm的系数；d(t r) i表示在第i观测历元测站接收机钟差的待求参数；
    则简化后的标准单点定位误差方程的线性化表达式，如公式(6)所示；
    Figure PCTCN2021123845-appb-100026
    公式(6)中，i＝1,2,…,epoch，epoch表示最大观测历元数；
    Figure PCTCN2021123845-appb-100027
    表示由卫星j在第i观测历元伪距观测值、卫星j在第i观测历元的钟差以及测站近似坐标与卫星j在第i观测历元 之间的几何距离计算出的常数项，其中，
    Figure PCTCN2021123845-appb-100028
    设定滑动窗口选择i个历元，每个历元最多能够观测到j颗卫星，则在这个滑动窗口下所有可视卫星的伪距观测值组成imax个方程，其中，imax＝i×j；
    定义标准单点定位误差方程的系数矩阵、伪距观测值的改正数向量矩阵、常数项矩阵以及未知参数矩阵分别为矩阵B、V、L、X，则它们的表达式分别如下：
    Figure PCTCN2021123845-appb-100029
    Figure PCTCN2021123845-appb-100030
    Figure PCTCN2021123845-appb-100031
    X＝[dX dY dZ (dt r) 1 … (dt r) i C 00 … C NmaxNmax S NmaxNmax] T；
    则标准单点定位误差方程的线性化表达式(6)用一个矩阵形式表达，如公式(7)所示：
    V＝BX-L  (7)
    未知参数X的最小二乘估值为：X＝(B TPB) -1B TPL  (8)
    公式(8)中，P为单位权矩阵；
    I.3.2.判断公式(7)中系数矩阵B是否为病态，经过判断，若系数矩阵B是病态，则执行步骤I.3.3；否则，系数矩阵B是良态，执行步骤I.3.4；
    I.3.3.首先利用截断奇异值分解的方法，为最小二乘谱修正迭代提供待求解的未知参数X的初值，然后基于最小二乘谱修正迭代计算未知参数X的值；
    设系数矩阵B∈R n×m，R n×m表示n行m列的实数阵；则系数矩阵B的奇异值分解为：
    B＝USV T  (9)
    公式(9)中，U∈n×n，V∈m×m，U、V均为正交矩阵，S∈n×m为一对角矩阵；
    系数矩阵B∈R n×m的截断奇异值矩阵B k定义为：
    Figure PCTCN2021123845-appb-100032
    公式(10)中，矩阵S中最小的(r-k)个非奇异值用零代替，即被截断，且k≤r；r表示系数矩阵B的秩，k表示矩阵S中保留的奇异值数量；
    u i表示矩阵U对应的向量，v i表示矩阵V对应的向量，σ i表示矩阵S中保留的奇异值；
    计算矩阵S中奇异值的均值，将该均值作为截断奇异值的阈值，其中，矩阵S中奇异值大于截断奇异值的阈值的保留，小于截断奇异值的阈值的进行零化处理；
    则公式(7)的截断奇异值解
    Figure PCTCN2021123845-appb-100033
    为：
    Figure PCTCN2021123845-appb-100034
    公式(11)中，
    Figure PCTCN2021123845-appb-100035
    根据最小二乘原理V TPV＝min，公式(8)写成以下形式：
    (B TPB)X＝(B TPL)  (12)
    将公式(12)的等式左端和右端同时加未知参数KX，并化简得到公式(13)：
    (B TPB+KI)X＝B TPL+KX  (13)
    由公式(13)整理得，最小二乘谱修正迭代公式为：
    X＝(B TPB+KI) -1(B TPL+KX)  (14)
    公式(14)中，K为任意实数；
    基于最小二乘谱修正迭代计算未知参数X的解，转到步骤I.3.5；
    I.3.4.利用最小二乘法求解未知参数X的解，转到步骤I.3.5；
    I.3.5.根据计算得到的未知参数X的解，得到未知参数X中包含的参数值，即测站的位置参数(dX,dY,dZ)、接收机钟差dt r以及球谐展开的系数C nm和S nm；
    将未知参数X中包含的位置参数(dX,dY,dZ)，分别加到测站的近似坐标(X 0,Y 0,Z 0)上，即标准单点定位求出的地心地固坐标系下的测站坐标(X,Y,Z)；
    I.3.6.将测站坐标由地心地固坐标系转换到站心坐标系，实现GNSS标准单点定位。
  2. 根据权利要求1所述基于球谐展开的GNSS标准单点定位方法，其特征在于，
    所述步骤I.3.2中，系数矩阵B是否为病态的判断过程如下：计算公式(7)的系数阵B的条件数，设置条件数的经验阈值；将条件数与经验阈值进行比较，经过比较判断：
    当条件数大于经验阈值时，则表明系数矩阵B为病态，否则，系数矩阵B为良态。
  3. 根据权利要求1所述基于球谐展开的GNSS标准单点定位方法，其特征在于，
    所述步骤I.3.3中，基于最小二乘谱修正迭代公式计算X的解的过程如下：
    基于最小二乘谱修正迭代解算X包括迭代①和迭代②两个迭代过程；其中，设定用于判断迭代①是否收敛的第一比较阈值，设定用于判断迭代②是否收敛的第二比较阈值；
    迭代①第一次迭代时，设定K的初始值为1，将公式(11)得到的截断奇异值解
    Figure PCTCN2021123845-appb-100036
    作为公式(14)第一次迭代时X的初值，并赋值给公式(14)等式右边的X；
    每次迭代过程中，将上一次迭代时利用公式(14)得到的未知参数X的值，作为本次迭代过程中X的初值，并赋值给公式(14)等式右边的X；
    基于公式(14)求解得到每次迭代过程中未知参数X的值；
    将每次迭代过程中求解得到的X的值，与每次迭代时X的初值进行差值运算；
    若经过差值运算，得到的差值大于第一比较阈值，则表示迭代不收敛，需要改正最小二乘谱修正迭代中的K值，即将K值递增1，继续上述迭代①过程；
    若经过差值运算得到的差值小于或等于第一比较阈值，则跳出迭代①，进入迭代②；
    迭代②将利用公式(14)求得的未知参数X的值与第二比较阈值进行比较；
    若利用公式(14)求解得到的未知参数X的值大于第二比较阈值，则表示迭代②未收敛，此时将位置参数(dX,dY,dZ)分别加到测站的近似坐标(X 0,Y 0,Z 0)上，更新测站的近似坐标，通过公式(5)对更新测站近似坐标后的标准单点定位误差方程进行线性化，通过公式(6)简化线性化后的标准单点定位误差方程，重新进入迭代①求解未知参数X的值；
    若利用公式(14)求解得到的未知参数X的值小于或等于第二比较阈值，则表示迭代②收敛，跳出迭代；此时公式(14)求得的未知参数X的值，即未知参数X的解。
  4. 根据权利要求1所述基于球谐展开的GNSS标准单点定位方法，其特征在于，
    所述步骤I.3.4中，基于最小二乘法求解未知参数X的过程如下：
    设定用于判断迭代是否收敛的第三比较阈值；每次迭代过程中，利用公式(8)求解出未知参数X的值，将本次迭代得到的未知参数X的值与第三比较阈值进行比较；
    若利用公式(8)求解得到的未知参数X的值大于第三比较阈值，则表示迭代未收敛；此时，将位置参数(dX,dY,dZ)分别加到测站的近似坐标(X 0,Y 0,Z 0)上，更新测站的近似坐标，执行公式(5)-公式(8)的计算过程，并重新求解未知参数X的值；
    若利用公式(8)求解得到的未知参数X的值小于或等于第三比较阈值，则表示迭代收敛，跳出迭代；此时，利用公式(8)得到的未知参数X的值，即待求解的未知参数X的解。
PCT/CN2021/123845 2021-03-17 2021-10-14 一种基于球谐展开的gnss单点定位方法 WO2022048694A1 (zh)

Priority Applications (2)

Application Number Priority Date Filing Date Title
ZA2022/03722A ZA202203722B (en) 2021-03-17 2022-03-31 Gnss standard point positioning method based on spherical harmonics
US17/730,567 US20220299652A1 (en) 2021-03-17 2022-04-27 Gnss standard point positioning method based on spherical harmonics

Applications Claiming Priority (2)

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

Related Child Applications (1)

Application Number Title Priority Date Filing Date
US17/730,567 Continuation US20220299652A1 (en) 2021-03-17 2022-04-27 Gnss standard point positioning method based on spherical harmonics

Publications (1)

Publication Number Publication Date
WO2022048694A1 true WO2022048694A1 (zh) 2022-03-10

Family

ID=76668228

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/CN2021/123845 WO2022048694A1 (zh) 2021-03-17 2021-10-14 一种基于球谐展开的gnss单点定位方法

Country Status (4)

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

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114935768A (zh) * 2022-07-13 2022-08-23 武汉大学 一种基于单基站的虚拟参考站的构建方法
CN115598673A (zh) * 2022-09-29 2023-01-13 同济大学（Cn） Igs gnss卫星钟差和轨道单天相邻产品边界处偏差计算方法
CN115616625A (zh) * 2022-10-08 2023-01-17 国家基础地理信息中心 一种gnss实时数据偏移方法及***
CN116029177A (zh) * 2023-03-14 2023-04-28 中国科学院空天信息创新研究院 一种融合球谐函数和精细网格的对流层延迟建模方法
CN116027459A (zh) * 2022-12-30 2023-04-28 桂林理工大学 基于数值天气预报数据的大气加权平均温度的计算方法
CN116609810A (zh) * 2023-05-19 2023-08-18 复旦大学 基于导航地基***的电离层四维电子密度动态预测方法

Families Citing this family (13)

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

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109613582A (zh) * 2018-12-17 2019-04-12 中国科学院国家授时中心 一种车载实时单频米级伪距定位方法
CN109709591A (zh) * 2018-12-07 2019-05-03 中国科学院光电研究院 一种面向智能终端的gnss高精度定位方法
CN113093242A (zh) * 2021-03-17 2021-07-09 山东科技大学 一种基于球谐展开的gnss单点定位方法

Family Cites Families (20)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102096084B (zh) * 2010-12-09 2012-12-26 东南大学 基于星间组合差分的精密单点定位方法
NL2009695C2 (en) * 2012-10-25 2014-05-06 Fugro N V Ppp-rtk method and system for gnss signal based position determination.
US9557419B2 (en) * 2012-12-18 2017-01-31 Trimble Inc. Methods for generating accuracy information on an ionosphere model for satellite navigation applications
CN103323888B (zh) * 2013-04-24 2015-06-17 东南大学 Gnss大气探测数据中对流层延迟误差的消除方法
NL2013472B1 (en) * 2014-09-15 2016-09-28 Fugro N V Integer Ambiguity-Fixed Precise Point Positioning method and system.
NL2013473B1 (en) * 2014-09-15 2016-09-28 Fugro N V Precise GNSS positioning system with improved ambiguity estimation.
EP3035080B1 (en) * 2014-12-16 2022-08-31 Trimble Inc. Navigation satellite system positioning involving the generation of correction information
EP3130943B1 (en) * 2015-08-14 2022-03-09 Trimble Inc. Navigation satellite system positioning involving the generation of tropospheric correction information
CN105182388B (zh) * 2015-10-10 2017-08-25 安徽理工大学 一种快速收敛的精密单点定位方法
CN106407560B (zh) * 2016-09-19 2019-03-19 武汉大学 表征大气各向异性的对流层映射函数模型的构建方法
CN107356947B (zh) * 2017-05-31 2019-06-18 中国科学院测量与地球物理研究所 基于单频导航卫星数据确定卫星差分伪距偏差的方法
CN107153209B (zh) * 2017-07-06 2019-07-30 武汉大学 一种短弧段低轨导航卫星实时精密定轨方法
CN107765275B (zh) * 2017-09-04 2019-12-27 深圳市时空导航科技有限公司 广域差分定位方法、装置、终端及计算机可读存储介质
CN107632313A (zh) * 2017-09-13 2018-01-26 航天恒星科技有限公司 基于相关性的卫星导航信号和sbas电文仿真方法
CN108845340A (zh) * 2018-06-01 2018-11-20 浙江亚特电器有限公司 基于gnss-rtk的定位方法
CN110275185B (zh) * 2019-07-11 2020-04-03 武汉大学 基于gnss和geo卫星的电离层投影函数建模方法
CN111551971B (zh) * 2020-05-14 2021-05-25 中国北方工业有限公司 一种支持异频gnss信号伪距差分定位的方法
CN111551975B (zh) * 2020-06-24 2023-05-09 辽宁工程技术大学 Bds/gps参考站低高度角卫星整周模糊度确定方法
CN112034500A (zh) * 2020-08-20 2020-12-04 上海华测导航技术股份有限公司 基于实时ppp模糊度固定技术的区域格网电离层建模方法
CN111812681B (zh) * 2020-08-24 2023-11-24 中国人民解放军海军工程大学 大气区域建模方法、装置、电子设备及存储介质

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109709591A (zh) * 2018-12-07 2019-05-03 中国科学院光电研究院 一种面向智能终端的gnss高精度定位方法
CN109613582A (zh) * 2018-12-17 2019-04-12 中国科学院国家授时中心 一种车载实时单频米级伪距定位方法
CN113093242A (zh) * 2021-03-17 2021-07-09 山东科技大学 一种基于球谐展开的gnss单点定位方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
FU ZHENG ET AL.: "Real-time single-frequency pseudorange positioning in China based on regional satellite clock and ionospheric models", GPS SOLUTIONS, 13 November 2019 (2019-11-13), XP036993686, ISSN: 1080-5370, DOI: 10.1007/s10291-019-0923-2 *
GUO JIN-YUN, LI GUO-WEI, KONG QIAO-LI, WANG SHU-YANG: "Modeling GPS multipath effect based on spherical cap harmonic analysis", TRANSACTIONS OF NONFERROUS METALS SOCIETY OF CHINA, ELSEVIER, AMSTERDAM, NL, vol. 24, no. 6, 1 June 2014 (2014-06-01), AMSTERDAM, NL , pages 1874 - 1879, XP055907557, ISSN: 1003-6326, DOI: 10.1016/S1003-6326(14)63266-0 *
LIU MINGLIANG, JINGPING CHEN, KELIANG DING, LONGHAO YU, YAJIE LIU: "Beidou / GPS Standard Single Point Positioning Ｒesult Analysis", JOURNAL OF BEIJING UNIVERSITY OF CIVIL ENGINEERING AND ARCHITECTURE, vol. 34, no. 4, 31 December 2018 (2018-12-31), pages 34 - 40, XP055907555, ISSN: 1004-6011, DOI: 10.19740/j.1004-6011.2018.04.06 *

Cited By (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114935768A (zh) * 2022-07-13 2022-08-23 武汉大学 一种基于单基站的虚拟参考站的构建方法
CN114935768B (zh) * 2022-07-13 2022-11-04 武汉大学 一种基于单基站的虚拟参考站的构建方法
CN115598673A (zh) * 2022-09-29 2023-01-13 同济大学（Cn） Igs gnss卫星钟差和轨道单天相邻产品边界处偏差计算方法
CN115598673B (zh) * 2022-09-29 2023-10-24 同济大学 Igs gnss卫星钟差和轨道单天相邻产品边界处偏差计算方法
CN115616625A (zh) * 2022-10-08 2023-01-17 国家基础地理信息中心 一种gnss实时数据偏移方法及***
CN116027459A (zh) * 2022-12-30 2023-04-28 桂林理工大学 基于数值天气预报数据的大气加权平均温度的计算方法
CN116027459B (zh) * 2022-12-30 2024-05-14 桂林理工大学 基于数值天气预报数据的大气加权平均温度的计算方法
CN116029177A (zh) * 2023-03-14 2023-04-28 中国科学院空天信息创新研究院 一种融合球谐函数和精细网格的对流层延迟建模方法
CN116609810A (zh) * 2023-05-19 2023-08-18 复旦大学 基于导航地基***的电离层四维电子密度动态预测方法
CN116609810B (zh) * 2023-05-19 2024-06-07 复旦大学 基于导航地基***的电离层四维电子密度动态预测方法

Also Published As

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

Similar Documents

Publication Publication Date Title
WO2022048694A1 (zh) 一种基于球谐展开的gnss单点定位方法
CN107102346B (zh) 一种基于北斗***的多天线测姿方法
CN106168672B (zh) 一种gnss多模单频rtk周跳探测方法及装置
CN109917356B (zh) 一种机载激光扫描***误差标定方法
US11867823B2 (en) System and method for determining GNSS positioning corrections
JP2010528320A (ja) リアルタイムキネマティック（ｒｔｋ）測位における距離依存性誤差の軽減
CN109901206B (zh) 一种基于低轨卫星无线电测距信号的单星定位与授时方法
CN106407560B (zh) 表征大气各向异性的对流层映射函数模型的构建方法
Kashcheyev et al. Estimation of higher-order ionospheric errors in GNSS positioning using a realistic 3-D electron density model
Yao et al. A new method to accelerate PPP convergence time by using a global zenith troposphere delay estimate model
CN113671505A (zh) 一种基于***几何误差补偿的合成孔径雷达立体定位方法
Zhou et al. Real-time orbit determination of Low Earth orbit satellite based on RINEX/DORIS 3.0 phase data and spaceborne GPS data
Wang et al. Performance evaluation of a real-time high-precision landslide displacement detection algorithm based on GNSS virtual reference station technology
Raquet et al. Use of a Covariance Analysis Technique for Predicting Performance of Regional‐Area Differential Code and Carrier‐Phase Networks
Li et al. Improve the ZY-3 height accuracy using ICESat/GLAS laser altimeter data
Wang et al. Instantaneous sub-meter level precise point positioning of low-cost smartphones
CN115980317B (zh) 基于修正相位的地基gnss-r数据土壤水分估算方法
CN104991265A (zh) 一种北斗卫星导航***用户统一性定位方法
Li et al. A grid-based ionospheric weighted method for PPP-RTK with diverse network scales and ionospheric activity levels
Gao et al. Ionosphere-weighted post-processing kinematic for airborne positioning with refined modeling of receiver phase biases and tropospheric zenith wet delays
CN115308781B (zh) 基于bdgim辅助的相位平滑伪距高精度时间传递方法
CN114910939B (zh) 短距离大高差rtk中对流层延迟实测气象改正方法
Zhang et al. BDS satellites and receivers DCB resolution
CN113465575A (zh) 一种基于对流层先验信息约束的高落差山区gnss高精度快速网解方法
Peng et al. Precise orbit determination for Jason-1 satellite using on-board GPS data with cm-level accuracy

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 21863756

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 21863756

Country of ref document: EP

Kind code of ref document: A1