CN111539109B - Real-time high-precision global multi-dimensional troposphere zenith delay grid model construction method - Google Patents

Real-time high-precision global multi-dimensional troposphere zenith delay grid model construction method Download PDF

Info

Publication number
CN111539109B
CN111539109B CN202010332446.XA CN202010332446A CN111539109B CN 111539109 B CN111539109 B CN 111539109B CN 202010332446 A CN202010332446 A CN 202010332446A CN 111539109 B CN111539109 B CN 111539109B
Authority
CN
China
Prior art keywords
ztd
window
global
grid
elevation
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202010332446.XA
Other languages
Chinese (zh)
Other versions
CN111539109A (en
Inventor
黄良珂
彭华
刘立龙
黎峻宇
康传利
谢劭峰
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Guilin University of Technology
Original Assignee
Guilin University of Technology
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 Guilin University of Technology filed Critical Guilin University of Technology
Priority to CN202010332446.XA priority Critical patent/CN111539109B/en
Publication of CN111539109A publication Critical patent/CN111539109A/en
Application granted granted Critical
Publication of CN111539109B publication Critical patent/CN111539109B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Position Fixing By Use Of Radio Waves (AREA)

Abstract

The invention belongs to the technical field of intersection of satellite navigation positioning system data processing and meteorology, and discloses a real-time high-precision global multi-dimensional troposphere zenith delay grid model construction method, wherein a sliding window algorithm is introduced to divide a global section into 5 degrees multiplied by 4 degrees, ERA5 data is utilized to construct a ZTD vertical section model considering the fine seasonal change of a ZTD elevation scaling factor in each window of the global, and then a5 degrees multiplied by 4 degrees global ZTD vertical section grid model is constructed. In the invention, multidimensional ZTD model parameters at the average elevation position of each window in the world are used as parameters of the geometric center of each window, and are stored in a grid point form, a new global grid is constructed according to a sliding window algorithm, and a real-time high-precision global multidimensional ZTD grid model suitable for any height is constructed. The real-time high-precision global multidimensional ZTD grid model established by the method well improves the calculation efficiency of the model and enhances the practicability of the model.

Description

Real-time high-precision global multi-dimensional troposphere zenith delay grid model construction method
Technical Field
The invention belongs to the technical field of Global Navigation Satellite System (GNSS) data processing and meteorology intersection, and particularly relates to a real-time high-precision global multi-dimensional troposphere zenith delay grid model construction method.
Background
Currently, tropospheric delay is a key factor affecting high-precision navigation positioning of spatial technologies such as Global Navigation Satellite System (GNSS), Very Long Baseline Interferometry (VLBI), and the like, and is also basic data of atmospheric scientific research. Tropospheric delay models are currently mainly used to provide tropospheric delay information. Traditional tropospheric delay models are mainly classified into two types, one is a model requiring meteorological parameters, and the other is a non-meteorological parameter model. Tropospheric delay models requiring actual measurement of meteorological parameters have poor real-time performance, thus limiting their application in real-time navigational positioning. The non-meteorological parameter troposphere delay model is widely applied to real-time navigation positioning because the meteorological parameters are not required to be measured actually and the use is convenient. The non-meteorological parameter troposphere delay correction models widely applied globally at present mainly include an EGNOS model, a UNB series model and the like. The models are global average tropospheric atmosphere delay models which only can reflect the global outline of the global tropospheric atmosphere space-time change, are difficult to reflect regional tropospheric atmosphere change characteristics, and have low space-time resolution.
Aiming at the defects, a plurality of scholars develop the construction research of the troposphere delay grid global model, such as a Tropsrid series model, a GPT series model and the like. In the TropGrid series model, the use precision of the TropGrid model is improved by 25 percent compared with that of the EGNOS model in the global scope. Tropigrid 2 is an improved version of the tropigrid model, which improves the temporal resolution of the model, but the model ignores the half-year-cycle variation in tropospheric delay. For a GPT series model, the troposphere delay information of the user position is acquired by providing meteorological parameters such as air temperature, air pressure, water vapor pressure and the like and combining an empirical troposphere delay model. Although the GPT series model is widely applied to space positioning technologies such as GNSS and VLBI, ECMWF monthly average profile data is adopted in modeling of the GPT series model, and daily period change of meteorological parameters is difficult to describe. In summary, the vertical resolution and temporal resolution of the model still need to be improved. In addition, with the continuous deepening of the research and application of the troposphere delay model abroad, domestic researchers also carry out corresponding research to construct various global troposphere delay models, such as an IGGtrop series model, a GZTD model, an ITG model and the like.
Through the above analysis, the problems and defects of the prior art are as follows: (1) the spatial resolution of the tropospheric vertical section model is not fine enough; (2) the model construction uses a single data source.
The difficulty in solving the above problems and defects is: constructing a ZTD vertical section grid model with high space-time resolution; and fusing multi-source data to construct a global multi-dimensional real-time ZTD grid model.
The significance of solving the problems and the defects is as follows: the method can provide high-precision real-time troposphere delay information at any height for satellite navigation positioning and space atmosphere detection, and promotes the development and application of GNSS meteorology.
Disclosure of Invention
Aiming at the problems in the prior art, the invention provides a real-time high-precision global multi-dimensional troposphere zenith delay grid model construction method.
The invention is realized in such a way that a real-time high-precision global multi-dimensional troposphere zenith delay grid model construction method comprises the following steps:
dividing the global into regular windows with the same size based on a sliding window algorithm, constructing a global ZTD vertical section grid model taking account of the fine seasonal variation of a ZTD elevation scaling factor by utilizing ERA5 data, fusing high-precision and high-space-time-resolution troposphere zenith delay ZTD information of a mesoscale weather forecast center, a global geodetic observation system atmospheric center, an IGS center and an urban CORS network platform, and constructing a global ZTD model taking account of the elevation, the latitude and the seasonal variation of each window in the global;
and taking the constructed parameters of the global window ZTD model as the parameters of the geometric center of each window, storing the parameters in a grid point form, constructing a global grid under a new resolution according to a sliding window algorithm, and constructing a real-time high-precision global multidimensional ZTD grid model.
The method specifically comprises the following steps:
acquiring the latest ERA5 data provided by an ECMWF center, grid ZTD data provided by a GGOS atmosphere center and the corresponding grid surface elevation, ZTD products provided by an IGS center and GNSS data provided by a CORS network as data sources, and performing corresponding data processing;
step two, dividing the global into a regular window of 5 degrees multiplied by 4 degrees based on a sliding window algorithm, and calculating ZTD information of each ERA5 grid point on different isostatic pressing layers in the global by utilizing the extracted global ERA5 grid point layered meteorological data with the annual planar resolution of 1 degree multiplied by 1 degree and corresponding surface meteorological data;
expressing a global ZTD vertical section function model of each window by using a negative index function, and estimating coefficients of the ZTD vertical section model of each window in the world, namely 5 coefficients of beta factors, by using ZTD section information of 30 ERA5 lattice points in each window with the resolution of 1 hour for many years according to each window in the world through a least square method;
step four, storing the obtained 5 coefficients of the beta factor of each window of the globe in a grid form with the plane resolution of 5 degrees multiplied by 4 degrees (longitude multiplied by latitude), and obtaining a ZTD vertical section grid model with the plane resolution of 5 degrees multiplied by 4 degrees (longitude multiplied by latitude) and taking account of the fine seasonal change of the ZTD elevation scaling factor;
step five, resampling different time resolution ZTD values which are precisely extracted from GNSS reference stations, namely IGS station 5min ZTD products and CO RS network GNSS station 1h ZTD data contained in each window of the world to obtain a ZTD value which has the same time resolution as GGOS grid ZTD, namely 6 h;
step six, performing elevation reduction on the GGOS grid ZTD product and the GNSS reference station ZTD data, and correcting the GGOS grid ZTD product and the GNSS reference station ZTD data to average elevation H in respective windows of the whole world0At least one of (1) and (b);
step seven, projecting multisource ZTD values with 6h resolution in each window of the world to the average elevation position of each window, and constructing a ZTD time sequence model of the internal view, elevation, latitude and fine seasonal change of each window of the world based on the fused multisource ZTD data;
and step eight, taking the multidimensional ZTD model parameters at the average height of each window in the world as the parameters of the geometric center of each window, storing the parameters in a grid point form, constructing a new global grid with the plane resolution of 5 degrees multiplied by 4 degrees (longitude multiplied by latitude) according to a sliding window algorithm, and constructing a real-time high-precision global multidimensional ZTD grid model suitable for any height.
Further, in the first step, the data source and the data processing method include:
the latest ERA5 data is the grid point hierarchical data and the surface data of the global ERA5 data with the multi-year horizontal resolution of 1 degree multiplied by 1 degree; each grid point is 37 layers of data; the layered data comprises air pressure, potential height, specific humidity and temperature data; the earth surface data comprise earth surface air pressure, earth surface temperature, earth surface specific humidity data and corresponding earth surface elevation;
aiming at GNSS reference station data provided by a CORS network, ZTD information of each GNSS reference station with resolution better than 1h is calculated by Bernese high-precision GNSS processing software.
Further, in step two, the rule window includes:
a total of 3240 windows; each window contains 30 ERA5 data grid dots of 1 ° × 1 ° planar resolution.
Further, in the second step, the ZTD information calculation formula includes:
Figure BDA0002465444590000041
Figure BDA0002465444590000042
Figure BDA0002465444590000043
wherein N represents the total refractive index of the atmosphere, P represents the atmospheric pressure (hPa), e represents the vapor pressure (hPa), Q represents the specific humidity, T represents the temperature, H represents the elevation, H _ low represents the bottommost height integrated by the ERA5 reanalysis data, H _ top represents the topmost height integrated by the ERA5 reanalysis data, k represents the total refractive index of the atmosphere, P represents the atmospheric pressure (hPa), e represents the vapor pressure (hPa), Q represents the specific humidity, T represents the temperature, H represents the elevation1＝77.604K/Pa、k264.79K/Pa and K3＝375463K2the/hPa is a constant coefficient; n is a radical ofiAnd Δ HiRespectively representing the total atmospheric refractive index and the atmospheric thickness of the ith layer, and n represents the integral atmospheric layer number; ZHDSaasERA5 representing Saastamoinen model calculations reanalyzed residual ZHD values at the top of the data layer;
the ZHDSaasThe calculation formula is as follows:
Figure BDA0002465444590000044
in the formula (I), the compound is shown in the specification,
Figure BDA0002465444590000045
is the latitude of the survey station, and the unit is rad; h is elevation in km.
Further, in step three, the estimating, by a least square method, the ZTD vertical section model coefficients for each window around the world includes:
the ZTD vertical profile function for each window is as follows:
ZTDHu＝ZTDHr*exp(-βi*(Hu-Hr))；
in the formula, Hu and Hr are expressed as target elevation and reference elevation, ZTDHuAnd ZTDHrThe ZTD values at the target elevation and the reference elevation, respectively, β is the ZTD elevation scaling factor, i denotes the number of the window (i ═ 1,2,3, …, 3240);
β for each window around the world, which can be expressed as:
Figure BDA0002465444590000051
in the formula, A0Is the annual mean value of ZTD elevation scaling factor (A)1，A2) Represents the annual cycle amplitude coefficient of the ZTD elevation scaling factor, (A)3，A4) Represents the half-year-cycle amplitude coefficient of the ZTD elevation scaling factor.
Further, in step six, the correcting the GGOS grid ZTD product and GNSS reference station ZTD data to average elevations within respective windows around the world comprises:
the average elevation of each window is the average elevation of the GGOS grid points in the window and the elevation of the GNSS reference station falling into the window;
the calculation formula of the average elevation of each window is as follows:
Figure BDA0002465444590000052
in the formula H0The average elevation of the window is represented,
Figure BDA0002465444590000053
representing the elevation of the ith GGOS grid point within the window,
Figure BDA0002465444590000054
representing the elevation, n, of the kth GNSS station within the window1Indicates the number of GGOS grid points in the window, n2Representing the number of GNSS reference stations within the window.
Further, in the seventh step, the ZTD time sequence model expression for the consideration of elevation, latitude and fine seasonal change in each window of the world is as follows:
Figure BDA0002465444590000055
Figure BDA0002465444590000056
Figure BDA0002465444590000057
in the formula (I), the compound is shown in the specification,
Figure BDA0002465444590000058
representing the mean elevation of the window H0A ZTD value of (a); ZTDHA ZTD value representing H at surface height; Δ ZTD represents the ZTD correction projected from the surface height to the average elevation of the window, and can be calculated by the constructed ZTD vertical section function in the window; a is0Is the ZTD annual mean at the window mean height;
Figure BDA0002465444590000059
representing the latitude; a is1Representing a latitude correction factor; a is2、a3And a4Respectively representing the amplitude of the ZTD year, half year and day periods at the window average height; phi is a1、φ2And phi3Phase of ZTD year, half year and day periods at the window average height are respectively represented; doy is the product date of the year, and hod is UTC, all multi-source data in the window are utilized, and the least square method is adopted to solve the ZTD time sequence model parameters at the average elevation in each window all over the world.
Further, the real-time high-precision global multi-dimensional troposphere zenith delay grid model construction method further comprises the following steps:
searching and acquiring a multidimensional ZTD model parameter of a grid point at the geometric center of a window closest to a user side according to the approximate position information of the user side; calculating the average height of the window where the user terminal is located according to the obtained multidimensional ZTD model parameters
Figure BDA0002465444590000061
Then, using window ZTD vertical section grid model
Figure BDA0002465444590000062
The values are reduced to ZTD values at the user end height.
Another objective of the present invention is to provide a real-time high-precision global multi-dimensional troposphere zenith delay grid model construction system, which comprises:
the different isobaric layer ZTD information acquisition module is used for acquiring the latest ERA5 data provided by the ECMWF center, the grid ZTD data provided by the GGOS atmosphere center and the corresponding grid surface elevation, the ZTD product provided by the IGS center and the GNSS data provided by the CORS network as data sources and carrying out corresponding data processing;
dividing the global into a regular window of 5 degrees multiplied by 4 degrees based on a sliding window algorithm, and calculating ZTD information of each ERA5 lattice point on different isostatic pressing layers in the global by utilizing extracted global ERA5 lattice point hierarchical meteorological data with annual planar resolution of 1 degree multiplied by 1 degree and corresponding surface meteorological data;
the coefficient acquisition module of the ZTD vertical profile function expresses a ZTD vertical profile function model of each global window by using a negative index function, and for each global window, the coefficients of the ZTD vertical profile function of each global window, namely 5 coefficients of a beta factor, are estimated by using ZTD profile information of 30 ERA5 lattice points in the window with the resolution of 1 hour for many years through a least square method;
the ZTD vertical section grid model acquisition module is used for storing the obtained 5 coefficients of the beta factor of each window in the world in a grid form with the plane resolution of 5 degrees multiplied by 4 degrees to obtain a ZTD vertical section grid model with the plane resolution of 5 degrees multiplied by 4 degrees and taking the fine seasonal change of the ZTD elevation scaling factor into consideration;
the ZTD value acquisition module is used for resampling ZTD values with different time resolutions, which are precisely extracted from GNSS reference stations, namely an IGS station 5minZTD product and a CORS network GNSS station 1h ZTD data, contained in each window of the world, and acquiring the ZTD value with the same time resolution as that of a GGOS grid ZTD, namely 6 h;
the ZTD elevation reduction module is used for carrying out elevation reduction on the GGOS grid ZTD product and the GNSS reference station ZTD data and correcting the GGOS grid ZTD product and the GNSS reference station ZTD data to the average elevation H in the respective windows of the whole world0At least one of (1) and (b);
the ZTD time sequence model building module is used for projecting multisource ZTD values with 6h resolution in each window of the world to the average elevation position of each window, and building the ZTD time sequence model for considering elevation, latitude and fine seasonal change in each window of the world based on the fused data;
the real-time high-precision global multidimensional ZTD grid model acquisition module takes multidimensional ZTD model parameters at the average height of each window in the world as parameters of the geometric center of each window, stores the parameters in a grid point form, constructs a new global grid with the plane resolution of 5 degrees multiplied by 4 degrees according to a sliding window algorithm, and constructs a real-time high-precision global multidimensional ZTD grid model suitable for any height.
Another object of the present invention is to provide a program storage medium for receiving user input, the stored computer program causing an electronic device to execute the method for constructing a real-time high-precision global multi-dimensional troposphere zenith delay grid model, comprising: dividing the global into regular windows with the same size based on a sliding window algorithm, constructing a global ZTD vertical section grid model taking account of the fine seasonal variation of a ZTD elevation scaling factor by utilizing ERA5 data, fusing high-precision and high-space-time-resolution troposphere zenith delay ZTD information of a mesoscale weather forecast center, a global geodetic observation system atmospheric center, an IGS center and an urban CORS network platform, and constructing a global ZTD model taking account of the elevation, the latitude and the seasonal variation of each window in the global;
and taking the constructed parameters of the global window ZTD model as the parameters of the geometric center of each window, storing the parameters in a grid point form, constructing a global grid under a new resolution according to a sliding window algorithm, and constructing a real-time high-precision global multidimensional ZTD grid model.
It is another object of the present invention to provide a computer program product stored on a computer readable medium, comprising a computer readable program for providing a user input interface to implement the real-time high-precision global multi-dimensional tropospheric zenith delay mesh model construction method when executed on an electronic device.
The invention also utilizes a global troposphere delay model and common troposphere products to carry out global precision inspection and regional applicability evaluation on the real-time high-precision global multidimensional ZTD refinement model.
By combining all the technical schemes, the invention has the advantages and positive effects that:
the real-time high-precision global multidimensional ZTD grid model established by the invention realizes the fusion modeling of multi-source data and solves the defect of low spatial resolution of the ZTD vertical section model. The real-time high-precision global multidimensional ZTD grid model established by the invention can provide high-precision ZTD information at any global position in real time and can also be used for real-time high-precision ZTD elevation reduction. Compared with the current GPT2w model with the resolution of 1 degree multiplied by 1 degree and the optimal nominal precision, the model parameters of the real-time high-precision global multidimensional ZTD grid model established by the invention are reduced by about 43 times. In addition, as the constructed global multidimensional ZTD grid model fuses GNSS station ZTD data, the regional enhancement of the model can be embodied when the model is used in a region with GNSS stations. In a word, the real-time high-precision global multidimensional ZTD grid model established by the method well improves the calculation efficiency of the model and the practicability of the model.
The invention provides a method for constructing a real-time high-precision global multidimensional ZTD grid new model fusing multi-source data on the basis of dividing a global section into regular windows with consistent sizes based on a sliding window algorithm, fusing troposphere zenith delay (ZTD) information with high precision and high space-time resolution provided by platforms such as an European mesoscale weather forecast center (ECMWF), a global geodetic observation system Atmosphere center (GGOS Atmosphere), an IGS center, a city CORS grid and the like, and constructing a vertical section function model taking account of the fine seasonal change of ZTD elevation scaling factors of each global window.
Drawings
Fig. 1 is a flowchart of a method for constructing a real-time high-precision global multi-dimensional troposphere zenith delay grid model according to an embodiment of the present invention.
Fig. 2 is a schematic diagram of a real-time high-precision global multi-dimensional troposphere zenith delay grid model construction method provided by the embodiment of the invention.
Fig. 3 is a schematic diagram of a sliding window algorithm provided in an embodiment of the present invention.
Detailed Description
In order to make the objects, technical solutions and advantages of the present invention more apparent, the present invention is further described in detail with reference to the following embodiments. It should be understood that the specific embodiments described herein are merely illustrative of the invention and are not intended to limit the invention.
Aiming at the problems in the prior art, the invention provides a real-time high-precision global multi-dimensional troposphere zenith delay grid model construction method, and the invention is described in detail below by combining the attached drawings.
As shown in fig. 1-2, a method for constructing a real-time high-precision global multi-dimensional troposphere zenith delay grid model according to an embodiment of the present invention includes:
s101, acquiring the latest ERA5 data provided by the ECMWF center, grid ZTD data provided by the GGOS atmosphere center and the grid surface elevation corresponding to the grid ZTD data, ZTD products provided by the IGS center and GNSS data provided by the CORS network as data sources, and performing corresponding data processing.
S102, dividing the global into a regular window of 5 degrees multiplied by 4 degrees based on a sliding window algorithm, and calculating ZTD information of each ERA5 grid point on different isobaric layers in the global by utilizing the extracted hierarchical meteorological data of the global ERA5 grid point with the annual plane resolution of 1 degree multiplied by 1 degree and the corresponding surface meteorological data.
S103, expressing a global ZTD vertical section function model of each window by using a negative index function, and estimating the ZTD vertical section model coefficients of each window, namely 5 coefficients of beta factors, of each window by using the ZTD section information of 30 ERA5 lattice points in each window with the resolution of 1 hour for many years according to each window in the world through a least square method.
And S104, storing the obtained 5 coefficients of each global window elevation scaling beta factor in a grid form with the plane resolution of 5 degrees multiplied by 4 degrees (longitude multiplied by latitude), and obtaining a ZTD vertical section grid model with the plane resolution of 5 degrees multiplied by 4 degrees (longitude multiplied by latitude) and considering the fine seasonal variation of the ZTD elevation scaling factor.
And S105, resampling different time resolution ZTD values obtained by precisely extracting 5min ZTD products of the GNSS reference station, namely the IGS station, and 1h ZTD data of the GNSS station of the CORS network in each window of the world, and obtaining the ZTD value with the same time resolution as that of the GGOS grid ZTD, namely 6 h.
S106, performing elevation reduction on the GGOS grid ZTD product and the GNSS reference station ZTD data, and correcting the GGOS grid ZTD product and the GNSS reference station ZTD data to average elevation H in respective windows of the whole world0To (3).
S107, projecting multisource ZTD values with 6h resolution in each window of the world to the average elevation of each window, and constructing a ZTD time sequence model of the global internal view, elevation, latitude and fine seasonal change based on the fused multisource ZTD data.
And S108, taking the multidimensional ZTD model parameters at the average height of each window in the world as the parameters of the geometric center of each window, storing the parameters in a grid point form, constructing a new global grid with the plane resolution of 5 degrees multiplied by 4 degrees (longitude multiplied by latitude) according to a sliding window algorithm, and constructing a real-time high-precision global multidimensional ZTD grid model suitable for any height.
In step S101, the data source and the data processing method provided in the embodiment of the present invention include:
the latest ERA5 data is the grid point hierarchical data and the surface data of the global ERA5 data with the multi-year horizontal resolution of 1 degree multiplied by 1 degree; each grid point is 37 layers of data; the layered data comprises air pressure, potential height, specific humidity and temperature data; the surface data comprises surface air pressure, surface temperature, surface specific humidity data and corresponding surface elevation.
Aiming at GNSS reference station data provided by a CORS network, ZTD information of each GNSS reference station with resolution better than 1h is calculated by Bernese high-precision GNSS processing software.
In step S102, 3240 rule windows provided by the embodiment of the present invention; each window contains 30 ERA5 data grid dots of 1 ° × 1 ° planar resolution.
In step S102, the ZTD information calculation formula provided in the embodiment of the present invention includes:
Figure BDA0002465444590000101
Figure BDA0002465444590000102
Figure BDA0002465444590000103
wherein N represents the total refractive index of the atmosphere, P represents the atmospheric pressure (hPa), e represents the vapor pressure (hPa), Q represents the specific humidity, T represents the temperature, H represents the elevation, H _ low represents the bottommost height integrated by the ERA5 reanalysis data, H _ top represents the topmost height integrated by the ERA5 reanalysis data, k represents the total refractive index of the atmosphere, P represents the atmospheric pressure (hPa), e represents the vapor pressure (hPa), Q represents the specific humidity, T represents the temperature, H represents the elevation1＝77.604K/Pa、k264.79K/Pa and K3＝375463K2the/hPa is a constant coefficient; n is a radical ofiAnd Δ HiRespectively representing the total atmospheric refractive index and the atmospheric thickness of the ith layer, and n represents the integral atmospheric layer number; ZHDSaasERA5, representing Saastamoinen model calculations, reanalyzed residual ZHD values at the top of the data layer.
The ZHDSaasThe calculation formula is as follows:
Figure BDA0002465444590000104
in the formula (I), the compound is shown in the specification,
Figure BDA0002465444590000111
is the latitude of the survey station, and the unit is rad; h is elevation in km.
In step S103, the estimating of the ZTD vertical section model coefficient of each global window by the least square method provided in the embodiment of the present invention includes:
the ZTD vertical profile function for each window is as follows:
ZTDHu＝ZTDHr*exp(-βi*(Hu-Hr))；
in the formula, Hu and Hr are expressed as target elevation and reference elevation, ZTDHuAnd ZTDHrThe ZTD values at the target elevation and the reference elevation, respectively, β is the ZTD elevation scaling factor, and i denotes the window number (i ═ 1,2,3, …, 3240).
β for each window around the world, which can be expressed as:
Figure BDA0002465444590000112
in the formula, A0Is the annual mean value of ZTD elevation scaling factor (A)1，A2) Represents the annual cycle amplitude coefficient of the ZTD elevation scaling factor, (A)3，A4) Represents the half-year-cycle amplitude coefficient of the ZTD elevation scaling factor.
In step S106, the correcting the GGOS grid ZTD product and the GNSS reference station ZTD data to the average elevation within each window in the world according to the embodiment of the present invention includes:
and the average elevation of each window is the average elevation of the GGOS grid points in the window and the elevation of the GNSS reference station falling into the window.
The calculation formula of the average elevation of each window is as follows:
Figure BDA0002465444590000113
in the formula H0The average elevation of the window is represented,
Figure BDA0002465444590000114
representing the elevation of the ith GGOS grid point within the window,
Figure BDA0002465444590000115
representing the elevation, n, of the kth GNSS station within the window1Indicates the number of GGOS grid points in the window, n2Representing the number of GNSS reference stations within the window.
In step S107, the ZTD time series model expression for the internal view of each window in the world, the elevation, the latitude, and the fine seasonal variation provided in the embodiment of the present invention is as follows:
Figure BDA0002465444590000116
Figure BDA0002465444590000121
Figure BDA0002465444590000122
in the formula (I), the compound is shown in the specification,
Figure BDA0002465444590000123
representing the mean elevation of the window H0A ZTD value of (a); ZTDHA ZTD value representing H at surface height; Δ ZTD represents the ZTD correction projected from the surface height to the average elevation of the window, and can be calculated by the constructed ZTD vertical section function in the window; a is0Is the ZTD annual mean at the window mean height;
Figure BDA0002465444590000124
representing the latitude; a is1Representing a latitude correction factor; a is2、a3And a4Respectively representing the amplitude of the ZTD year, half year and day periods at the window average height; phi is a1、φ2And phi3Phase of ZTD year, half year and day periods at the window average height are respectively represented; doy is annual product date, hod is UTC, all multi-source data in the window are utilized, and the least square method is adopted to solve the ZTD time sequence model parameter at the average height in each window in the world.
The use method of the real-time high-precision global multi-dimensional troposphere zenith delay grid model provided by the embodiment of the invention comprises the following steps:
searching and acquiring a multidimensional ZTD model parameter of a grid point at the geometric center of a window closest to a user according to the approximate position information of the user; calculating the ZTD at the average height of the window where the user is located according to the obtained multidimensional ZTD model parametersH0Using window ZTD vertical section grid model to convert ZTD into three-dimensional grid modelH0The values are reduced to ZTD values at the user height.
The method for testing the real-time high-precision global multi-dimensional troposphere zenith delay grid model provided by the embodiment of the invention comprises the following steps:
and carrying out global precision inspection and regional applicability evaluation on the real-time high-precision global multidimensional ZTD refinement model by utilizing a global troposphere delay model and a common troposphere product.
The invention provides a real-time high-precision global multi-dimensional troposphere zenith delay grid model construction system, which comprises:
and the different isobaric layer ZTD information acquisition module is used for acquiring the latest ERA5 data provided by the ECMWF center, the grid ZTD data provided by the GGOS atmosphere center and the corresponding grid surface elevation, the ZTD product provided by the IGS center and the GNSS data provided by the CORS network as data sources and carrying out corresponding data processing.
Based on a sliding window algorithm, the global is divided into a regular window of 5 degrees multiplied by 4 degrees, and the extracted global ERA5 grid dot hierarchical meteorological data with the annual planar resolution of 1 degree multiplied by 1 degree and the corresponding surface meteorological data are utilized to calculate the ZTD information of each worldwide ERA5 grid dot on different equal-pressure layers.
The coefficient acquisition module of the ZTD vertical profile function expresses a ZTD vertical profile function model of each global window by using a negative index function, and for each global window, the coefficients of the ZTD vertical profile function of each global window, namely 5 coefficients of a beta factor, are estimated by using ZTD profile information of 30 ERA5 lattice points in the window with the resolution of 1 hour for many years through a least square method.
And the ZTD vertical section grid model acquisition module is used for storing the obtained 5 coefficients of the beta factor of each window in the world in a grid form with the plane resolution of 5 degrees multiplied by 4 degrees to obtain the ZTD vertical section grid model with the plane resolution of 5 degrees multiplied by 4 degrees and taking the fine seasonal change of the ZTD elevation scaling factor into consideration.
And the ZTD value acquisition module is used for resampling ZTD values with different time resolutions after precisely extracting GNSS reference stations, namely an IGS station 5minZTD product and CORS network GNSS station 1h ZTD data contained in each window of the world, and acquiring the ZTD value with the same time resolution as that of the GGOS grid ZTD, namely 6 h.
The ZTD elevation reduction module is used for carrying out elevation reduction on the GGOS grid ZTD product and the GNSS reference station ZTD data and correcting the GGOS grid ZTD product and the GNSS reference station ZTD data to the average elevation H in the respective windows of the whole world0To (3).
And the ZTD time sequence model building module is used for projecting multisource ZTD values with 6h resolution in each window of the world to the average elevation position of each window, and building the ZTD time sequence model for considering elevation, latitude and fine seasonal change in each window of the world based on the fused data.
The real-time high-precision global multidimensional ZTD grid model acquisition module takes multidimensional ZTD model parameters at the average height of each window in the world as parameters of the geometric center of each window, stores the parameters in a grid point form, constructs a new global grid with the plane resolution of 5 degrees multiplied by 4 degrees according to a sliding window algorithm, and constructs a real-time high-precision global multidimensional ZTD grid model suitable for any height.
The technical solution of the present invention is further illustrated by the following specific examples.
Example 1:
(1) data source introduction and processing. The model construction of the invention relates to data sources which are respectively as follows: the latest ERA5 data provided by ECMWF center, grid ZTD data provided by GGOS atmosphere center and its corresponding grid surface elevation, ZTD product provided by IGS center and GNSS data provided by CORS grid. Wherein the ERA5 reanalyzed data has a horizontal resolution of 0.25 degree, a temporal resolution of 1 hour, and a vertical resolution of 37 equal-pressure layers; the horizontal resolution of the GGOS large grid ZTD product is 2.5 degrees multiplied by 2 degrees (longitude multiplied by latitude), and the time resolution is 6 h; the IGS center provides a ZTD product of thousands of GNSS reference stations in the world, the time resolution is 5min, and the precision is better than 5 mm; aiming at GNSS reference station data provided by a CORS network, calculating ZTD information of each GNSS reference station with a resolution ratio superior to 1h by using Bernese high-precision GNSS processing software; the method needs to download and extract the hierarchical data of the grid points (each grid point is 37 layers of data) of the global ERA5 data with the horizontal resolution of 1 degree multiplied by 1 degree and the earth surface data for many years, wherein the hierarchical data comprises the air pressure, the potential height, the specific humidity and the temperature data, and the earth surface data comprises the earth surface air pressure, the earth surface temperature, the earth surface specific humidity data and the corresponding earth surface elevation; the global ZTD vertical section grid model is constructed by using ZTD values calculated by ERA5 data, and on the basis, GGOS large grid ZTD products and GNSS reference station ZTD data are fused to construct a global multidimensional ZTD grid model.
(2) The principle of ZTD calculation. The ZTD value of each grid point in the world at the surface elevation and the ZTD values of different equal-pressure layers can be obtained by integral calculation, and the integral formula for calculating the ZTD is as follows:
Figure BDA0002465444590000141
Figure BDA0002465444590000142
Figure BDA0002465444590000143
wherein N represents the total refractive index of the atmosphere, P represents the atmospheric pressure (hPa), e represents the vapor pressure (hPa), Q represents the specific humidity,t represents temperature, H represents elevation, H _ low represents the bottommost height calculated by the ERA5 reanalysis data integration, H _ top represents the topmost height calculated by the ERA5 reanalysis data integration, k1＝77.604K/Pa、k264.79K/Pa and K3＝375463K2The values of/hPa are constant coefficients. Since the ERA5 reanalyzes the top layer of the data and residual atmosphere, in order to improve the computing precision of ZTD, a Saastamoinen model is adopted to compute the residual tropospheric delay value above the layer top, and the residual tropospheric delay value is added to the integration result of each final layer of each lattice point of ERA 5. The integral calculation is convenient, the formula (1) needs to be discretized, and can be rewritten as follows:
Figure BDA0002465444590000144
wherein N isiAnd Δ HiRespectively representing the total atmospheric refractive index and the atmospheric thickness of the ith layer, and n represents the integral atmospheric layer number; ZHDSaasERA5, representing the Saastamoinen model calculation, reanalyzes the residual ZHD values at the top of the data layer, as calculated by the following equation:
Figure BDA0002465444590000151
in the formula (I), the compound is shown in the specification,
Figure BDA0002465444590000152
the latitude (unit: rad) and the elevation (unit: km) of the measuring station are shown.
(3) A sliding window algorithm is introduced. The determination of the size of the sliding window needs to take the principles of the integer of the global window subdivision number, the continuity of the window, the resolvability of model parameters in the window and the like into consideration. Based on the above principle, the present invention exemplifies the sliding window algorithm with the area range of 5 ° × 4 ° (longitude × latitude) as a sliding window size, and the flow is shown in fig. 3. The specific process is as follows: firstly, utilizing the first window N at the upper left corner of the grid1Inner (each black box represents a sliding window size) multi-source data obtains troposphere related model parameters in the window, and the troposphere related model parameters are used as the windowMouth N1The result of the center grid point (black dot in box); then the window is moved to the east direction of the latitude by 2 lattice points, and a new window N is solved2Inner troposphere-related model parameters, as window N2Calculating the troposphere related model parameters in all windows at the latitude by analogy according to the result of the central grid point; then the window is moved to the next latitude (two grid points are moved downwards), the troposphere related model parameters in all windows of the latitude are solved in the same method, and the like is repeated until the troposphere related model parameters in all windows of the world are solved. Finally, the troposphere-related model parameters in all the windows of the world are obtained and taken as the results of the central grid points of the respective windows, and finally, a new global grid is constructed from the central grid points of all the windows of the world, as shown by the black points and the dotted lines in fig. 3.
(3) And (4) performing integral calculation on information of each window ZTD vertical section of the world. The invention divides the world into regular windows (3240 windows) of 5 degrees multiplied by 4 degrees based on a sliding window algorithm, each window comprises 30 ERA5 data grid points with 1 degree multiplied by 1 degree plane resolution, and the ZTD information of each ERA5 grid point on different isobaric layers in the world is calculated according to the integration of the formulas (1) - (5) by utilizing the extracted global ERA5 grid point layered meteorological data (air pressure, potential height, specific humidity and temperature) with the annual plane resolution of 1 degree multiplied by 1 degree and the corresponding surface meteorological data (air pressure, surface elevation, specific humidity and temperature).
(4) And constructing a ZTD vertical section model. Expressing a global ZTD vertical profile function model of each window by using a negative exponential function, and estimating a ZTD vertical profile model coefficient of each window by using ZTD profile information of 30 ERA5 grid points in each window with the resolution of 1 hour for many years according to a least square method for each window in the world, wherein the ZTD vertical profile function of each window is as follows:
ZTDHu＝ZTDHr*exp(-βi*(Hu-Hr)) (6)
in the formula, Hu and Hr are expressed as target elevation and reference elevation, ZTDHuAnd ZTDHrThe ZTD values at the target elevation and the reference elevation, respectively, β is the ZTD elevation scaling factor, i denotes the window number (i ═ 1,2, 3)…, 3240). The 5-year time sequence of the beta factor is subjected to spectrum detection analysis by adopting Fast Fourier Transform (FFT), and the beta factor is mainly represented by annual cycle and semiannual cycle changes. To this end, β for each window around the world, which can be expressed as:
Figure BDA0002465444590000161
in the formula, A0Is the annual mean value of ZTD elevation scaling factor (A)1，A2) Represents the annual cycle amplitude coefficient of the ZTD elevation scaling factor, (A)3，A4) Represents the half-year-cycle amplitude coefficient of the ZTD elevation scaling factor. The 5 coefficients of the global beta factor of each window are stored in a grid form with the plane resolution of 5 degrees multiplied by 4 degrees (longitude multiplied by latitude), and a ZTD vertical section grid model which takes account of the fine seasonal variation of the ZTD elevation scaling factor and has the plane resolution of 5 degrees multiplied by 4 degrees (longitude multiplied by latitude) is obtained.
(5) The ZTD vertical profile mesh model was used. The use of the new model is very simple, with the following steps: 1) searching a lattice point closest to the ZTD vertical profile lattice model according to the plane position (longitude and latitude) of the user, and acquiring model parameters (5 coefficients of a beta factor) of the lattice point; 2) based on the elevation information of the user, the ZTD at the reference elevation can be finally determined according to equation (6)HrZTD of values to elevation normalization to target heightHuThe value is obtained.
(6) And (5) testing the precision of the ZTD vertical section grid model. And (3) verifying the accuracy of the global ZTD vertical section grid model by combining ZTD information of different equal-pressure layers of the global grid point which does not participate in the modeling ERA5 data integral calculation and ZTD information of different sections of the global probe station data integral calculation.
(7) The multi-source ZTD data is resampled in time. In order to realize the multi-source data joint modeling, firstly, the time resolutions of different data sources are unified, so that the time resolution ZTD values precisely extracted by GNSS reference stations (IGS station 5minZTD product and CORS net GNSS station 1h ZTD data) contained in various windows around the world are resampled to obtain the ZTD value with the same time resolution (6h) as that of the GGOS grid ZTD.
(8) All-purposeAnd calculating the average elevation of each window of the ball. In order to realize the fusion modeling of multi-source data, the elevation reduction of the GGOS grid ZTD product and the GNSS reference station ZTD data is needed to be carried out, and the elevation reduction is corrected to the average elevation H in each window of the whole world0And the average elevation of each window is the average elevation of the GGOS grid point in the window and the elevation of the GNSS reference station falling into the window, and the calculation formula is as follows:
Figure BDA0002465444590000171
in the formula H0The average elevation of the window is represented,
Figure BDA0002465444590000172
representing the elevation of the ith GGOS grid point within the window,
Figure BDA0002465444590000173
representing the elevation, n, of the kth GNSS station within the window1Indicates the number of GGOS grid points in the window, n2Representing the number of GNSS reference stations within the window.
(9) And fusing multi-source data to construct a global ZTD time sequence model with high time resolution of each window. And projecting the multisource ZTD value with 6h resolution in each window of the world to the average elevation of each window to realize multisource data fusion modeling. The ZTD fine seasonal variation detection is the basis for constructing a high-time-resolution ZTD time sequence model, and research shows that the ZTD shows remarkable annual cycle and half-annual cycle variation in the world. Therefore, the inventor conducts deep research on the day-cycle variation characteristic of global ZTD by using GGOS grid ZTD, and shows that the ZTD also has obvious day-cycle variation, and particularly in middle and low latitude areas, the day-cycle variation of the ZTD is more obvious. In addition, the ERA5 data is used for preliminary study on the spatial change of the ZTD, and the study shows that the ZTD is influenced by the elevation direction much more than the plane direction, and the ZTD has smaller change on different longitudes on the same latitude band at the same height. Based on the analysis, the invention constructs a ZTD time sequence model of the internal view of each window of the world, the elevation, the latitude and the fine seasonal variation, and the expression formula is shown as the formula (9) and the formula (10):
Figure BDA0002465444590000174
Figure BDA0002465444590000175
in the formula (I), the compound is shown in the specification,
Figure BDA0002465444590000176
representing the mean elevation of the window H0A ZTD value of (a); ZTDHA ZTD value representing H at surface height; the delta ZTD represents the ZTD correction value projected from the earth surface height to the window average elevation, and can be calculated through a constructed ZTD vertical section model in the window; a is0Is the ZTD annual mean at the window mean height;
Figure BDA0002465444590000177
representing the latitude; a is1Representing a latitude correction factor; a is2、a3And a4Respectively representing the amplitude of the ZTD year, half year and day periods at the window average height; phi is a1、φ2And phi3Phase of ZTD year, half year and day periods at the window average height are respectively represented; doy is year day, hod is UTC. And solving the ZTD time sequence model parameter at the average height in each window in the world by using all multi-source data in the window and adopting a least square method.
(10) And constructing a global multi-dimensional ZTD grid model with high precision and high space-time resolution. Global fine expression of multi-dimensional ZTD model parameters (ZTD time sequence model parameters, ZTD vertical profile model parameters and global average elevation of each window) of each window is the key for constructing a high space-time resolution global multi-dimensional ZTD refinement model. In the invention, multidimensional ZTD model parameters at the average height of each window in the world are used as parameters of the geometric center of each window, and are stored in a grid point form, and a new global grid with the plane resolution of 5 degrees multiplied by 4 degrees (longitude multiplied by latitude) is finally constructed according to a sliding window algorithm, and a real-time high-precision global multidimensional ZTD grid model applicable to any height is finally constructed.
(11) And (3) using and testing a global multi-dimensional ZTD refinement model with high space-time resolution. And searching and acquiring a multidimensional ZTD model parameter of a grid point at the geometric center of the window closest to the user according to the approximate position information of the user. According to the obtained multidimensional ZTD model parameters, firstly, calculating the average height of the window where the user is positioned
Figure BDA0002465444590000181
Then using window ZTD vertical section grid model
Figure BDA0002465444590000182
The values are reduced to ZTD values at the user height.
Meanwhile, a common global troposphere delay model (GPT2w, UNB3m, EGNOS and other models) and common troposphere products (products such as ZTD provided by an IGS center, ZTD provided by a GGOS atmosphere center, ZTD calculated by CORS network reference station data, ERA5 data, radio detection empty station data and the like) are used for carrying out global precision inspection and regional applicability evaluation on the real-time high-precision global multidimensional ZTD refinement model constructed by the method, and the effectiveness, reliability and regional enhancement of the model are verified. The general technical route is shown in fig. 2.
Multisource data are resampled, elevation reduction is carried out on the multisource ZTD data based on the constructed ZTD vertical section model, and the multisource ZTD data are reduced to the average elevation of each window; establishing a ZTD time sequence model of the internal consideration of each window, the elevation, the latitude and the fine seasonal change of the whole world based on the multisource ZTD data of each window after the elevation is reduced; and taking the multidimensional ZTD model parameters at the average elevation position of each window in the world as the parameters of the geometric center of each window, storing the parameters in a grid point form, constructing a new global 5-degree multiplied by 4-degree grid according to a sliding window algorithm, and constructing a real-time high-precision global multidimensional ZTD grid model suitable for any height. The real-time high-precision global multidimensional ZTD grid model established by the method well improves the calculation efficiency of the model and enhances the practicability of the model.
Through the above description of the embodiments, those skilled in the art will clearly understand that the present invention may be implemented by software plus a necessary hardware platform, and may also be implemented by hardware entirely. With this understanding in mind, all or part of the technical solutions of the present invention that contribute to the background can be embodied in the form of a software product, which can be stored in a storage medium, such as a ROM/RAM, a magnetic disk, an optical disk, etc., and includes instructions for causing a computer device (which can be a personal computer, a server, or a network device, etc.) to execute the methods according to the embodiments or some parts of the embodiments of the present invention.
The above description is only for the purpose of illustrating the present invention and the appended claims are not to be construed as limiting the scope of the invention, which is intended to cover all modifications, equivalents and improvements that are within the spirit and scope of the invention as defined by the appended claims.

Claims (8)

1. A real-time high-precision global multi-dimensional troposphere zenith delay grid model construction method is characterized by comprising the following steps:
dividing the global into regular windows with the same size based on a sliding window algorithm, constructing a ZTD vertical grid model taking account of the fine seasonal variation of a ZTD elevation scaling factor, fusing the high-precision and high-space-time-resolution troposphere zenith delay ZTD information of a mesoscale weather forecast center, a global geodetic observation system atmospheric center, an IGS center and a city CORS network platform, and constructing a ZTD model taking account of the elevation, the latitude and the seasonal variation of each window in the global;
taking the constructed parameters of each global window ZTD model as the parameters of the geometric center of each window, storing the parameters in a grid point form, constructing a global grid under a new resolution according to a sliding window algorithm, and constructing a real-time high-precision global multidimensional ZTD grid model;
the real-time high-precision global multi-dimensional troposphere zenith delay grid model construction method further comprises the following steps:
acquiring the latest ERA5 data provided by an ECMWF center, grid ZTD data provided by a GGOS atmosphere center and the corresponding grid surface elevation, ZTD products provided by an IGS center and GNSS data provided by a CORS network as data sources, and performing corresponding data processing;
step two, dividing the global into a regular window of 5 degrees multiplied by 4 degrees based on a sliding window algorithm, and calculating ZTD information of each ERA5 grid point on different isostatic pressing layers in the global by utilizing the extracted global ERA5 grid point layered meteorological data with the annual planar resolution of 1 degree multiplied by 1 degree and corresponding surface meteorological data;
expressing a global ZTD vertical section function model of each window by using a negative index function, and estimating coefficients of the global ZTD vertical section model of each window by using ZTD section information of 30 ERA5 lattice points in each window with the resolution of 1 hour for many years according to each window in the world through a least square method, wherein the coefficients are 5 coefficients of beta factors;
step four, storing the obtained 5 coefficients of the beta factor of each window in the world in a grid form with the plane resolution of 5 degrees multiplied by 4 degrees to obtain a ZTD vertical section grid model with the plane resolution of 5 degrees multiplied by 4 degrees and taking the fine seasonal change of the ZTD elevation scaling factor into consideration;
step five, resampling different time resolution ZTD values which are precisely extracted from GNSS reference stations, namely IGS station 5min ZTD products and C ORS network GNSS station 1h ZTD data contained in each window of the world to obtain a ZTD value which has the same time resolution as the GGOS grid ZTD, namely 6 h;
step six, performing elevation reduction on the GGOS grid ZTD product and the GNSS reference station ZTD data, and correcting the GGOS grid ZTD product and the GNSS reference station ZTD data to average elevation H in respective windows of the whole world0At least one of (1) and (b);
step seven, based on the established global ZTD vertical section grid model, projecting multisource ZTD values with 6h resolution in each window of the world to the average elevation position of each window, and constructing a ZTD time sequence model for the internal consideration and elevation, latitude and fine seasonal change of each window of the world based on the fused multisource ZTD data;
and step eight, taking the multi-dimensional ZTD model parameters at the average elevation position of each window in the world as the parameters of the geometric center of each window, storing the parameters in a grid point mode, constructing a new global grid with the plane resolution of 5 degrees multiplied by 4 degrees according to a sliding window algorithm, and constructing a real-time high-precision global multi-dimensional ZTD grid model suitable for any height.
2. The method for constructing the real-time high-precision global multi-dimensional troposphere zenith delay grid model according to claim 1, wherein in the first step, the data source and the data processing method comprise:
the latest ERA5 data is the grid point hierarchical data and the surface data of the global ERA5 data with the multi-year horizontal resolution of 1 degree multiplied by 1 degree; each grid point is 37 layers of data; the layered data comprises air pressure, potential height, specific humidity and temperature data; the earth surface data comprise earth surface air pressure, earth surface temperature, earth surface specific humidity data and corresponding earth surface elevation;
aiming at GNSS reference station data provided by a CORS network, calculating ZTD information of each GNSS reference station with a resolution ratio superior to 1h by using Bernese high-precision GNSS processing software;
in the second step, 3240 rule windows are provided; each window contains 30 ERA5 data grid dots of 1 ° × 1 ° planar resolution;
in the second step, the ZTD information calculation formula includes:
Figure FDA0002997038000000021
Figure FDA0002997038000000022
Figure FDA0002997038000000023
wherein N represents the total refractive index of the atmosphere, P represents the atmospheric pressure, and the unit is hPa; e represents the water vapour pressure in hPa; q represents specific humidityT represents temperature, H represents elevation, H _ low represents the bottommost height calculated by the ERA5 reanalysis data integration, H _ top represents the topmost height calculated by the ERA5 reanalysis data integration, k1＝77.604K/Pa、k264.79K/Pa and K3＝375463K2the/hPa is a constant coefficient; n is a radical ofiAnd Δ HiRespectively representing the total atmospheric refractive index and the atmospheric thickness of the ith layer, and n represents the integral atmospheric layer number; ZHDSaasERA5 representing Saastamoinen model calculations reanalyzed residual ZHD values at the top of the data layer;
the ZHDSaasThe calculation formula is as follows:
Figure FDA0002997038000000031
in the formula (I), the compound is shown in the specification,
Figure FDA0002997038000000032
is the latitude of the survey station, and the unit is rad; h is elevation in km.
3. The method for constructing the real-time high-precision global multi-dimensional troposphere zenith delay grid model of claim 1, wherein in step three, the estimating the coefficients of the ZTD vertical section model of each global window by the least square method comprises:
the ZTD vertical profile function for each window is as follows:
ZTDHu＝ZTDHr*exp(-βi*(Hu-Hr))；
in the formula, Hu and Hr are expressed as target elevation and reference elevation, ZTDHuAnd ZTDHrThe ZTD values at the target elevation and the reference elevation, respectively, β is the ZTD elevation scaling factor, i represents the number of the window, where i is 1,2, 3.
β for each window around the world, expressed as:
Figure FDA0002997038000000033
in the formula, A0Is the annual mean value of the ZTD elevation scaling factor, A1，A2Respectively representing the annual cycle amplitude coefficient, A, of the ZTD elevation scaling factor3，A4Respectively represent the half-year-cycle amplitude coefficients of the ZTD elevation scaling factors.
4. The method for constructing a real-time high-precision global multi-dimensional troposphere zenith delay mesh model of claim 1 wherein, in step six, the correcting the GGOS mesh ZTD products and GNSS reference station ZTD data to average elevations within respective windows around the world comprises:
the average elevation of each window is the average elevation of the GGOS grid points in the window and the elevation of the GNSS reference station falling into the window;
the calculation formula of the average elevation of each window is as follows:
Figure FDA0002997038000000041
in the formula H0The average elevation of the window is represented,
Figure FDA0002997038000000042
representing the elevation of the ith GGOS grid point within the window,
Figure FDA0002997038000000043
representing the elevation, n, of the kth GNSS station within the window1Indicates the number of GGOS grid points in the window, n2Representing the number of GNSS reference stations within the window.
5. The method for constructing the real-time high-precision global multi-dimensional troposphere zenith delay grid model of claim 1, wherein in step seven, the ZTD time sequence model expression of the global window considering elevation, latitude and fine seasonal change is as follows:
Figure FDA0002997038000000044
Figure FDA0002997038000000045
in the formula (I), the compound is shown in the specification,
Figure FDA0002997038000000046
representing the mean elevation of the window H0A ZTD value of (a); ZTDHA ZTD value representing H at surface height; the delta ZTD represents the ZTD correction value projected from the earth surface height to the window average elevation, and can be calculated through a constructed ZTD vertical section model in the window; a is0Is the ZTD annual mean at the window mean height;
Figure FDA0002997038000000049
representing the latitude; a is1Representing a latitude correction factor; a is2、a3And a4Respectively representing the amplitude of the ZTD year, half year and day periods at the window average height; phi is a1、φ2And phi3Phase of ZTD year, half year and day periods at the window average height are respectively represented; doy is annual product date, hod is UTC, all multi-source data in the window are utilized, and the least square method is adopted to solve the ZTD time sequence model parameter at the average height in each window in the world.
6. The method for constructing the real-time high-precision global multi-dimensional troposphere zenith delay grid model according to claim 1, wherein the method for constructing the real-time high-precision global multi-dimensional troposphere zenith delay grid model further comprises:
searching and acquiring a multidimensional ZTD model parameter of a grid point at the geometric center of a window closest to a user side according to the approximate position information of the user side; calculating the average height of the window where the user terminal is located according to the obtained multidimensional ZTD model parameters
Figure FDA0002997038000000047
Then, using window ZTD vertical section grid model
Figure FDA0002997038000000048
The values are reduced to ZTD values at the user end height.
7. A real-time high-precision global multi-dimensional troposphere zenith delay grid model construction system of the real-time high-precision global multi-dimensional troposphere zenith delay grid model construction method according to any one of claims 1 to 6, wherein the real-time high-precision global multi-dimensional troposphere zenith delay grid model construction system comprises:
the different isobaric layer ZTD information acquisition module is used for acquiring the latest ERA5 data provided by the ECMWF center, the grid ZTD data provided by the GGOS atmosphere center and the corresponding grid surface elevation, the ZTD product provided by the IGS center and the GNSS data provided by the CORS network as data sources and carrying out corresponding data processing;
dividing the global into a regular window of 5 degrees multiplied by 4 degrees based on a sliding window algorithm, and calculating ZTD information of each ERA5 lattice point on different isostatic pressing layers in the global by utilizing extracted global ERA5 lattice point hierarchical meteorological data with annual planar resolution of 1 degree multiplied by 1 degree and corresponding surface meteorological data;
the coefficient acquisition module of the ZTD vertical profile function expresses the ZTD vertical profile model of each global window by using a negative index function, and for each global window, the coefficients of the ZTD vertical profile function of each global window, namely 5 coefficients of a beta factor, are estimated by using ZTD profile information of 30 ERA5 lattice points in the window with the resolution of 1 hour for many years through a least square method;
the ZTD vertical section grid model acquisition module is used for storing the obtained 5 coefficients of the beta factor of each window in the world in a grid form with the plane resolution of 5 degrees multiplied by 4 degrees to obtain a ZTD vertical section grid model with the plane resolution of 5 degrees multiplied by 4 degrees and taking the fine seasonal change of the ZTD elevation scaling factor into consideration;
the ZTD value acquisition module is used for resampling ZTD values with different time resolutions after precisely extracting GNSS reference stations, namely IGS station 5min ZTD products and CORS network GNSS station 1h ZTD data contained in each window of the world, and acquiring the ZTD value with the same time resolution as GGOS grid ZTD, namely 6 h;
the ZTD elevation reduction module is used for carrying out elevation reduction on the GGOS grid ZTD product and the GNSS reference station ZTD data and correcting the GGOS grid ZTD product and the GNSS reference station ZTD data to the average elevation H in the respective windows of the whole world0At least one of (1) and (b);
the ZTD time sequence model building module is used for projecting multisource ZTD values with 6h resolution in each window of the world to the average elevation position of each window, and building a ZTD time sequence model for the internal view of each window of the world and the changes of elevation, latitude and fine seasons based on the fused multisource ZTD data;
the real-time high-precision global multidimensional ZTD grid model acquisition module takes multidimensional ZTD model parameters at the average elevation position of each window in the world as parameters of the geometric center of each window, stores the parameters in a grid point mode, constructs a new global grid with the plane resolution of 5 degrees multiplied by 4 degrees according to a sliding window algorithm, and constructs a real-time high-precision global multidimensional ZTD grid model suitable for any height.
8. A program storage medium for receiving user input, the stored computer program causing an electronic device to execute the method for constructing a real-time high-precision global multi-dimensional troposphere zenith delay grid model of any one of claims 1 to 7, comprising: dividing the global into regular windows with the same size based on a sliding window algorithm, constructing a global ZTD vertical section grid model taking account of the fine seasonal variation of a ZTD elevation scaling factor by utilizing ERA5 data, fusing high-precision and high-space-time-resolution troposphere zenith delay ZTD information of a mesoscale weather forecast center, a global geodetic observation system atmospheric center, an IGS center and an urban CORS network platform, and constructing a global ZTD model taking account of the elevation, the latitude and the seasonal variation of each window in the global;
and taking the constructed parameters of the global window ZTD model as the parameters of the geometric center of each window, storing the parameters in a grid point form, constructing a global grid under a new resolution according to a sliding window algorithm, and constructing a real-time high-precision global multidimensional ZTD grid model.
CN202010332446.XA 2020-04-24 2020-04-24 Real-time high-precision global multi-dimensional troposphere zenith delay grid model construction method Active CN111539109B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010332446.XA CN111539109B (en) 2020-04-24 2020-04-24 Real-time high-precision global multi-dimensional troposphere zenith delay grid model construction method

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010332446.XA CN111539109B (en) 2020-04-24 2020-04-24 Real-time high-precision global multi-dimensional troposphere zenith delay grid model construction method

Publications (2)

Publication Number Publication Date
CN111539109A CN111539109A (en) 2020-08-14
CN111539109B true CN111539109B (en) 2021-05-21

Family

ID=71975430

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010332446.XA Active CN111539109B (en) 2020-04-24 2020-04-24 Real-time high-precision global multi-dimensional troposphere zenith delay grid model construction method

Country Status (1)

Country Link
CN (1) CN111539109B (en)

Families Citing this family (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112034490B (en) * 2020-10-10 2022-03-22 山东科技大学 NWP inversion troposphere delay improvement method
CN112560314B (en) * 2020-12-18 2023-08-22 南京信息工程大学 Method for improving accuracy of GPT2w model in calculating meteorological parameters
CN113639893B (en) * 2021-06-29 2022-09-30 东南大学 Near-earth weighted average temperature information acquisition method based on multiple meteorological factors
CN113960635B (en) * 2021-10-25 2024-04-30 山东科技大学 Troposphere delay correction method considering daily variation
NL2032589B1 (en) * 2022-07-25 2023-06-28 Naval University Eng Pla Method for constructing zenith troposphere delay database based on era-5
CN116029177B (en) * 2023-03-14 2023-06-09 中国科学院空天信息创新研究院 Troposphere delay modeling method integrating spherical harmonic function and fine grid

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104777488A (en) * 2015-03-13 2015-07-15 中国科学院上海天文台 Modeling method and device for zenith tropospheric delay as well as measuring method and device

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103323888B (en) * 2013-04-24 2015-06-17 东南大学 Method for eliminating delay errors of troposphere of GNSS atmospheric probing data
CN104965207B (en) * 2015-05-19 2018-02-09 同济大学 A kind of acquisition methods of zone convection layer zenith delay

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104777488A (en) * 2015-03-13 2015-07-15 中国科学院上海天文台 Modeling method and device for zenith tropospheric delay as well as measuring method and device

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
A grid-based tropospheric product for China using a GNSS network;Hongxing Zhang等;《Journal of Geodesy》;20171129;第92卷;第765-777页 *
Development of an improved empirical model for slant delays in the troposphere (GPT2w);Johannes Bohm等;《GPS Solut》;20140826;第19卷;第433-441页 *
SHAtrop：基于陆态网GNSS数据的中国大陆区域ZTD模型;陈俊平等;《武汉大学学报•信息科学版》;20191130;第44卷(第11期);第1588-1595页 *
地基GNSS对流层天顶延迟改正模型与方法研究;黄良珂;《中国优秀硕士学位论文全文数据库 基础科学辑》;20150115(第01期);第1-85页 *

Also Published As

Publication number Publication date
CN111539109A (en) 2020-08-14

Similar Documents

Publication Publication Date Title
CN111539109B (en) Real-time high-precision global multi-dimensional troposphere zenith delay grid model construction method
Xu et al. Urban morphology detection and computation for urban climate research
Yu et al. Evaluation of a high-resolution historical simulation over China: climatology and extremes
CN103323888B (en) Method for eliminating delay errors of troposphere of GNSS atmospheric probing data
Huang et al. A global grid model for the correction of the vertical zenith total delay based on a sliding window algorithm
CN106324620A (en) Tropospheric zenith delay method based not on real-time measurement of surface meteorological data
Zhao et al. High-precision ZTD model of altitude-related correction
Kim et al. Statistical downscaling for daily precipitation in Korea using combined PRISM, RCM, and quantile mapping: Part 1, methodology and evaluation in historical simulation
Gegout et al. Adaptive mapping functions to the azimuthal anisotropy of the neutral atmosphere
CN108154193A (en) A kind of long-term sequence precipitation data NO emissions reduction method
Kuhlmann et al. A novel gridding algorithm to create regional trace gas maps from satellite observations
Maksyutov et al. A high-resolution inverse modelling technique for estimating surface CO 2 fluxes based on the NIES-TM–FLEXPART coupled transport model and its adjoint
CN111126466B (en) Multi-source PWV data fusion method
Ohunakin et al. Solar radiation variability in Nigeria based on multiyear RegCM3 simulations
CN110310370B (en) Method for point-plane fusion of GPS (Global positioning System) and SRTM (short Range TM)
CN111538943B (en) Novel high-space-time resolution global ZTD vertical section grid model construction method
CN113960635B (en) Troposphere delay correction method considering daily variation
CN116609859A (en) Weather disaster high-resolution regional mode forecasting system and method
KR102002593B1 (en) Method and apparatus for analyzing harmful gas diffusion in a specific space
CN113639893B (en) Near-earth weighted average temperature information acquisition method based on multiple meteorological factors
Liu et al. On the study of influences of different factors on the rapid tropospheric tomography
Huang et al. A global grid model for the estimation of zenith tropospheric delay considering the variations at different altitudes
Ma et al. Establishment of regional tropospheric delay model in Australia
Liu et al. A Two-step projected iterative algorithm for tropospheric water vapor tomography
Rahman et al. Simulation of Rainfall over Bangladesh Using Regional Climate Model (RegCM4. 7)

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
EE01 Entry into force of recordation of patent licensing contract

Application publication date: 20200814

Assignee: Guangxi Taihua Information Technology Co.,Ltd.

Assignor: GUILIN University OF TECHNOLOGY

Contract record no.: X2022450000082

Denomination of invention: Real time and high-precision global multi-dimensional tropospheric zenith delay grid model construction method

Granted publication date: 20210521

License type: Common License

Record date: 20221118

Application publication date: 20200814

Assignee: Guangxi Guigong surveying and mapping Geographic Information Technology Co.,Ltd.

Assignor: GUILIN University OF TECHNOLOGY

Contract record no.: X2022450000077

Denomination of invention: Real time and high-precision global multi-dimensional tropospheric zenith delay grid model construction method

Granted publication date: 20210521

License type: Common License

Record date: 20221118

EE01 Entry into force of recordation of patent licensing contract