GNSS (Global Navigation Satellite System)-based offshore troposphere atmospheric parameter profile inversion method

Abstract

The invention relates to the field of marine meteorological monitoring, and discloses a GNSS-based offshore troposphere atmospheric parameter profile inversion method, which comprises the following steps of: obtaining actual observation values of total troposphere zenith delay, dry delay and wet delay; obtaining troposphere atmospheric parameter profile statistical characteristics through clustering analysis; data points related to temperature, water vapor pressure and air pressure are constructed by utilizing troposphere atmospheric parameter profile statistical characteristics, and theoretical values of zenith total delay, dry delay and wet delay corresponding to the data points are calculated through forward modeling; and selecting one or more of total troposphere zenith delay, dry delay and wet delay, utilizing an optimization algorithm to invert a troposphere atmospheric parameter profile, performing iteration to obtain an optimal data point, performing data quality control on the temperature, the water vapor pressure, the air pressure and the atmospheric refractive index corresponding to the optimal data point, and outputting the temperature, the water vapor pressure, the air pressure and the atmospheric refractive index. The method disclosed by the invention is low in equipment installation and maintenance cost, is suitable for an offshore observation platform, is reliable in inversion result, and can provide support for offshore troposphere atmosphere research.

Description

Technical Field

[0001] The present invention relates to the field of marine meteorological monitoring, and in particular to a method for inverting marine tropospheric atmospheric parameter profiles based on GNSS. Background Art

[0002] The atmospheric parameters of the marine troposphere include tropospheric temperature, water vapor pressure, air pressure and atmospheric refractive index. The vertical distribution characteristics of tropospheric temperature, water vapor pressure, air pressure and other parameters provide key input parameters for the data assimilation system, and at the same time play an important role in improving the forecast accuracy of extreme weather events such as typhoons and thunderstorms that affect production and life at sea. In the field of military applications, the atmospheric refractive index of the troposphere directly affects the bending characteristics of the electromagnetic wave propagation path, and thus determines the formation and parameter characteristics of the atmospheric waveguide. The evaporation waveguide or surface waveguide phenomenon formed under specific meteorological conditions will cause abnormal propagation of radar waves, which will cause detection blind spots for shipborne radars and can also be used for over-the-horizon communications. Therefore, the accurate inversion of the tropospheric atmospheric refractive index profile plays an important role in accurately identifying the waveguide type and key parameters.

[0003] There are two types of methods for obtaining the atmospheric parameter profile of the troposphere: contact and non-contact. The contact method usually uses the method of launching sounding balloons and sounding rockets, and the non-contact method includes microwave radiometers, lidar, etc. There are certain shortcomings in the application of the above detection methods to offshore observation platforms. The time resolution of launching sounding balloons and sounding rockets is limited, and they need to be launched manually. The operability on offshore platforms is low, and continuous observation at sea cannot be achieved. Microwave radiometers can achieve higher time resolution, but due to limitations such as high requirements for the stability of the installation platform, high power consumption, and high installation and maintenance costs, they are difficult to apply to sea-based observation platforms. Detection by optical methods such as lidar is easily affected by the high humidity and high salt environment at sea and is not corrosion-resistant. Although the autonomous sounding observation technology using drones and unmanned ships has developed rapidly in recent years, in actual use, unmanned sounding systems in offshore environments need to solve the problem of replenishing sounding sensors. Summary of the invention

[0004] In order to solve the above technical problems, the present invention provides a method for inverting the atmospheric parameter profile of the offshore troposphere based on GNSS, so as to achieve the purpose of low equipment installation and maintenance cost, being suitable for offshore observation platforms, and having reliable inversion results, so as to provide support for the research on the offshore troposphere atmosphere.

[0005] To achieve the above object, the technical solution of the present invention is as follows: A method for inverting atmospheric parameter profiles of the marine troposphere based on GNSS comprises the following steps: Step 1: Use GNSS real-time single-point positioning to estimate the total tropospheric zenith delay of the measuring station, use meteorological element sensors to measure the meteorological elements of the measuring station, calculate the tropospheric zenith dry delay based on the meteorological elements, separate the tropospheric zenith dry delay from the tropospheric zenith total delay, and then obtain the tropospheric zenith wet delay, thereby obtaining the actual observation values ​​of the tropospheric zenith total delay, dry delay and wet delay; Step 2, cluster analysis is performed using historical sounding data from sounding stations near the station or reanalysis data at the station location to obtain statistical characteristics of regional tropospheric atmospheric parameter profiles; Step 3, using the statistical characteristics of the tropospheric atmospheric parameter profile to construct data points on temperature, water vapor pressure and air pressure, and using the ray tracing method or numerical integration method combined with the atmospheric refractive index to forward calculate the theoretical values ​​of the zenith total delay, dry delay and wet delay corresponding to each data point; Step 4, select one or more of the total zenith delay, dry delay and wet delay of the troposphere, use the optimization algorithm to invert the tropospheric atmospheric parameter profile, select the data point set with the smallest error between the theoretical value and the actual observation value for iteration, and finally obtain the optimal data point. The temperature, water vapor pressure, air pressure and atmospheric refractive index corresponding to the optimal data point are subjected to data quality control and output as the final inversion result.

[0006] In the above scheme, for stations with co-location observation conditions of ceilometers and microwave radiometers, or conditions for receiving output data of numerical prediction models, the output data of ceilometers, microwave radiometers or numerical prediction models are used as constraints of the tropospheric atmospheric parameter profile inversion algorithm; among them, the cloud height data measured by the ceilometer will be used to constrain the relative humidity to be close to saturation near the cloud layer height, with a negative gradient above the cloud layer and a positive gradient below the cloud layer; the temperature, air pressure and water vapor pressure profiles output by the microwave radiometer or numerical prediction model are used as the initial values ​​of the data constructed in the subsequent inversion algorithm; If the observation station does not have the conditions for co-location observation of ceilometer and microwave radiometer, and does not have the conditions for receiving output data of numerical prediction model, no constraints will be set.

[0007] In the above scheme, in step 1, the method for obtaining the actual observed values ​​of the total tropospheric zenith delay, dry delay and wet delay is as follows: (1) Using GNSS receivers to obtain dual-frequency raw observations, combined with Beidou satellite-based broadcast PPP-B2b corrections, the ionospheric delay error is eliminated through ionospheric-free combination, and the actual observation value of the total tropospheric zenith delay is extracted. ; (2) Calculate the tropospheric zenith dry delay using meteorological elements: ; in, is the actual observed value of the tropospheric zenith dry delay, is the geographical latitude of the site, is the air pressure at the measuring station, is the station elevation; (3) Calculate the tropospheric zenith wet delay using the tropospheric zenith total delay and the tropospheric zenith dry delay: ; in, is the actual observed value of the tropospheric zenith wet delay.

[0008] In the above scheme, the clustering algorithm is a self-organizing map clustering algorithm. The specific method of using the self-organizing map clustering algorithm to obtain the statistical characteristics of the tropospheric atmospheric parameter profile is as follows: Select the historical sounding data of the sounding station near the station or the reanalysis data of the station location to construct the historical profile data set; pre-process the temperature, relative humidity and refractive index profiles in the data set: first, detect and remove outliers in the original profile data, then calculate the range of meteorological parameter data in each altitude layer, the first-order change rate and second-order change rate of the single profile of the parameter with altitude for each vertical layer data, normalize the data, and unify the parameters and their change rates of different altitude layers; The normalized profile data set is grouped by month, and the data set is subjected to nonlinear dimensionality reduction and topological mapping. The SOM self-organizing mapping network is trained and iteratively adjusted through competitive learning and neighborhood function. For each input data, its Euclidean distance with each neuron node is calculated, and data with similar characteristics are clustered to the same neuron node. After the iteration termination condition is met, the training is completed, and each data is assigned to the cluster represented by the most similar neuron to obtain the final clustering result. For each cluster, its statistical characteristics are calculated. First, the cluster center is extracted, and then the value range of meteorological parameters is hierarchically statistically analyzed along the height dimension. Further, each meteorological element is extracted and its first-order change rate and second-order change rate are calculated to obtain the intra-class statistical distribution characteristics, and finally the statistical characteristics of the tropospheric atmospheric parameter profile are obtained.

[0009] In the above scheme, in step 3, the specific method of constructing data points about temperature, water vapor pressure and air pressure by using the statistical characteristics of the tropospheric atmospheric parameter profile is as follows: Assume that the constructed data points are ,in, Indicates temperature, represents the water vapor pressure, Indicates air pressure, Indicates data points, represents the temperature profile, represents the water vapor pressure profile, represents the pressure profile; For temperature profile , is the total number of sections, express The first layer of data Layer; Assume that the minimum value of a single layer of the section is , the maximum value is , based on the profile characteristics, the first-order rate of change of the profile between the extreme values ​​is , whose value range is limited to the minimum value With the maximum value The second-order rate of change is , whose value range is limited to the minimum value With the maximum value between; For the first layer of random numbers , and its generation formula is: ; Here, Representation interval uniformly distributed random variables on , thus ensuring Values ​​in the interval Inside; Subsequent random numbers , for the first-order rate of change limit, by defining the first-order rate of change of each layer , the deformation is: ; in, is the height difference between the upper and lower layers of data; For the second-order rate limit without considering the height variable, it is defined as , the deformation is: ; Take the intersection of the first-order and second-order rate of change constraints for the random numbers of the above two formulas, and finally get the one that satisfies all the constraints. The value range is: ; If the intersection is empty, take the random number generated for the first-order rate of change limit, that is, , Representation interval Uniformly distributed random variables on ; According to the above process, a temperature profile that meets the requirements is generated based on the given single-layer maximum value and the statistical characteristics of the random change rate of the profile. The random number sequence of Sequence and pressure profile sequence, composition as a data point.

[0010] In the above scheme, in step 3, the method of forward modeling using ray tracing is as follows: First, the position of the GNSS antenna at the measuring station is taken as the starting position, and the position vector , ray tracing is performed towards the zenith; the vertical gradient of the atmospheric refractive index is calculated according to the ray tracing step length Interpolate to obtain the atmospheric refractive index corresponding to each calculation step ; From Layer to Layer iterative calculation until the calculation step size The accumulated height reaches the upper limit of the inversion profile; Secondly, the delay of signal propagation is calculated based on the layered atmosphere model. For each segment on the ray path, the refractive index of the segment is used to calculate the delay of signal propagation. The propagation delay of the signal in this section is calculated by the path length, and then the delays of each section are accumulated to obtain the theoretical value of the zenith total delay of the upper limit of the inversion profile; above the upper limit of the profile, the water vapor is zero, and the zenith total delay is represented by the zenith dry delay through the Saastamoinen model; the forward calculation formula for the theoretical value of the total tropospheric propagation delay is as follows: ; in, is the theoretical value of the total delay at the zenith of the troposphere, is the refractive index of each layer of atmosphere, For the ray The path length in the layer, is the zenith dry delay above the upper limit height of the profile; ; in, is the air pressure at the upper limit height of the section, is the latitude at the upper limit of the profile, is the upper limit height of the section; Similarly, the theoretical value of the tropospheric zenith dry delay is calculated as: ; in, is the theoretical value of the tropospheric zenith dry delay, is the dry refractive index of each layer of atmosphere; Theoretical value of tropospheric zenith wet delay The difference between the theoretical value of the total delay at the zenith of the troposphere and the theoretical value of the dry delay at the zenith of the troposphere is obtained: .

[0011] In the above scheme, in step 3, the method of forward modeling using the numerical integration method is as follows: 1) Assuming that the propagation path of GNSS signals in the troposphere is a straight line, the propagation path of GNSS signals in the troposphere is discretized and divided into several segments. Based on the discretization of the propagation path, the penetration points of each layer are calculated, that is, the points where the GNSS signal propagation path intersects vertically with a specific height layer. The position information of the penetration points is expressed as: ; in, For the The position of the layer's penetration point, is the latitude coordinate of the crossing point, is the longitude coordinate of the crossing point, is the penetration point height; 2) Based on the constructed data points, calculate the atmospheric refractive index profile corresponding to each data point: ; in, Indicates Data point The temperature of the layer, Indicates Data point The water vapor pressure of the layer, Indicates Data point The air pressure of the layer, For the Data point The atmospheric refractive index of the layer is Total Layer, namely: ; 3) The atmospheric refractive index is accumulated along the signal propagation path by the path integration method to obtain the theoretical value of the total delay at the zenith of the troposphere The calculation formula is: ; Theoretical value of tropospheric zenith dry delay The calculation formula is: ; in, is the atmospheric dry refractive index, is the zenith dry delay above the upper limit of the profile, From the Saastamoinen model: ; in, Height of penetration point The air pressure at Height of penetration point The latitude of is the penetration point height; Theoretical value of tropospheric zenith wet delay The theoretical value of the total delay at the zenith of the troposphere Theoretical value of tropospheric zenith dry delay The difference is: .

[0012] In the above scheme, the and The calculation formula is as follows: ; in, is the latitude of the GNSS receiver, is the longitude of the GNSS receiver, For the The satellite azimuth corresponding to the layer-through-layer point, The calculation formula is as follows: ; in, For the The satellite elevation angle corresponding to the layer-through-layer point, is the radius of the Earth, is the GNSS receiver altitude, For the Height of layer penetration point.

[0013] In the above scheme, in step 4, the optimization algorithm includes a single-objective optimization method or a multi-objective optimization method; wherein the single-objective optimization adopts a genetic algorithm or a particle swarm optimization algorithm, and the optimization goal is to minimize the root mean square error between the theoretical value and the actual observed value of one of the tropospheric zenith total delay, the tropospheric zenith dry delay and the tropospheric zenith wet delay; The multi-objective optimization adopts non-dominated sorting genetic algorithm or multi-objective particle swarm algorithm, and its optimization goal is to minimize the root mean square error between the theoretical value and the actual observation value of at least two of the tropospheric zenith total delay, tropospheric zenith dry delay and tropospheric zenith wet delay, which can be expressed as follows: ; in, is the multi-objective function vector, is the root mean square error between the actual observed value and the theoretical value of the zenith total delay, is the root mean square error between the actual observed value and the theoretical value of the zenith interference delay, is the root mean square error between the actual observed value and the calculated value of the zenith wet delay, represents the actual observed value of the total delay, Represents the theoretical value of total delay, represents the actual observed value of the dry delay, represents the theoretical value of dry delay, represents the actual observed value of wet delay; Indicates the theoretical wet delay value.

[0014] In a further technical solution, in step 4, when setting the constraint conditions, the method of inverting the tropospheric atmospheric parameter profile using the optimization algorithm is as follows: Using Constraint Sets To represent the set of external constraints: Constraints are set using cloud height data measured by ceilometers: profile data is divided into Layers, each layer height is , defines the cloud height measured by the ceilometer Corresponding section layer number , get the layer number closest to the cloud height; according to the existing temperature and air pressure Get saturated water vapor pressure , then the constraint condition for: ; in, Indicates Data point The water vapor pressure of the layer, Indicates Data point The saturated vapor pressure of the layer, is the layer corresponding to the cloud height, Indicates the number of layers corresponding to the cloud height The water vapor pressure, Cloud height corresponding to the number of layers The saturated water vapor pressure; Using the atmospheric parameter profile output by microwave radiometer or numerical prediction model, select the external output profile that can be obtained in actual situation to set the constraint conditions; for the constructed data point profile With the output profile , set the absolute value of the difference between the temperature and water vapor pressure data of each layer of the two to be no greater than a certain minimum value is a constraint condition, then the constraint condition for: ; in, Indicates that the constructed data point is The temperature of the layer, Indicates the external output profile in The temperature of the layer, Indicates that the constructed data point is The water vapor pressure of the layer, Indicates the external output profile in The water vapor pressure of the layer, Indicates that the constructed data point is The air pressure of the layer, Indicates the external output profile in The air pressure of the layer, , and It is The error thresholds of temperature, air pressure and water vapor pressure allowed by the layer; The final objective function is expressed as: ; The point set closest to the actual observation value is selected from the inverted data point set. Through multiple rounds of iterations and gradual optimization, the optimal data points are obtained. The temperature, water vapor pressure, air pressure and atmospheric refractive index corresponding to the optimal data points are subjected to data quality control and then output as the inversion result.

[0015] Through the above technical solution, the GNSS-based marine tropospheric atmospheric parameter profile inversion method provided by the present invention has the following beneficial effects: 1. The present invention uses GNSS signals and meteorological sensors to coordinate observations, combines Beidou satellite-based broadcasting of PPP-B2b corrections and dual-frequency ionospheric error elimination technology, and uses the Saastamoinen model to obtain high-precision actual observation values ​​of tropospheric zenith total delay, dry delay, and wet delay; 2. The present invention adopts the self-organizing map (SOM) clustering algorithm to construct a statistical characteristic model of the tropospheric atmospheric parameter profile, construct initial values ​​for subsequent inversion, thereby improving the inversion accuracy; 3. The present invention uses ray tracing method or numerical integration method for forward calculation, taking into account the signal path bending effect and calculation efficiency, and obtains the theoretical values ​​of the total delay, dry delay and wet delay at the zenith of the troposphere, providing accurate theoretical data for setting the objective function of the subsequent inversion; 4. The present invention uses an optimization algorithm combined with dynamic external constraints for inversion, considers using a multi-objective function as a target, and sets external constraints for the inversion algorithm according to specific circumstances, further improving the accuracy of the inversion result, thereby obtaining the optimal solution. In summary, the tropospheric atmospheric parameter profile inversion method based on GNSS signals provided by the present invention can achieve tropospheric atmospheric parameter profile inversion that is independent of historical sounding data, has the unique advantage of being suitable for offshore observation platforms that lack sounding data, has low equipment installation and maintenance costs, and can provide support for offshore tropospheric atmospheric research. BRIEF DESCRIPTION OF THE DRAWINGS

[0016] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art are briefly introduced below.

[0017] Figure 1 A schematic flow chart of a method for inverting an atmospheric parameter profile of the marine troposphere based on GNSS disclosed in an embodiment of the present invention; Figure 2 This is a schematic diagram of the installation of the measuring station equipment involved in the present invention. DETAILED DESCRIPTION

[0018] The technical solutions in the embodiments of the present invention will be described clearly and completely below in conjunction with the accompanying drawings in the embodiments of the present invention.

[0019] The present invention provides a method for inverting the atmospheric parameter profile of the marine troposphere based on GNSS, such as Figure 1 As shown, the following steps are included: Step 1: Use GNSS real-time single-point positioning to estimate the total tropospheric zenith delay of the measuring station location, use meteorological element sensors to measure the meteorological elements of the measuring station location, calculate the tropospheric zenith dry delay based on the meteorological elements, separate the tropospheric zenith dry delay from the tropospheric zenith total delay, and then obtain the tropospheric zenith wet delay, thereby obtaining the actual observation values ​​of the tropospheric zenith total delay, dry delay and wet delay.

[0020] The measuring station equipment of the present invention is installed as follows Figure 2 As shown, the main equipment involved includes: meteorological element sensors, GNSS antennas, GNSS receivers, data processing units, and external equipment such as cloud meters. Among them, meteorological element sensors and GNSS antennas generally need to be installed at the same height so as not to block each other, and their relative positions can be appropriately adjusted according to the specific conditions of the installation platform. The present invention provides a method for inverting the atmospheric parameter profile of the troposphere at sea based on GNSS, and its input data is collected from meteorological element sensors and GNSS receivers, and the algorithm software of the method runs on the data processing unit. For measuring stations that already have meteorological element sensors and GNSS receivers, it is only necessary to add a data processing unit and corresponding algorithm software to more conveniently realize the inversion function of the atmospheric parameter profile of the troposphere, and the equipment installation and maintenance costs are low.

[0021] The method to obtain the actual observed values ​​of the total tropospheric zenith delay, dry delay and wet delay is as follows: (1) Using GNSS receivers to obtain dual-frequency raw observations, combined with Beidou satellite-based broadcast PPP-B2b corrections, the ionospheric delay error is eliminated through ionospheric-free combination, and the actual observation value of the total tropospheric zenith delay is extracted. ; Air pressure at GNSS antenna location ,temperature and relative humidity Meteorological elements such as GNSS are measured by meteorological element sensors; GNSS raw observation pseudorange , carrier phase , Doppler shift and carrier-to-noise ratio Obtained by a GNSS receiver, the correction is obtained using a GNSS receiver or an external calculation unit.

[0022] The signal processing board is used to process the signal transmitted by the GNSS antenna, convert the signal into a digital signal, and complete the signal capture by searching for satellite signals through correlation operations in the digital domain. The navigation data in the signal is demodulated, and the binary data sequence containing parameters such as satellite orbit, clock, ionosphere model, etc. is extracted according to the reverse operation of the modulation method and transmitted to the data processing module.

[0023] Use high-precision GNSS data processing software RTKLIB for real-time precise point positioning. Receive and parse observation data and navigation messages from GNSS receivers, and perform quality control and preprocessing on pseudorange and carrier phase information in the observation data. Through precise point positioning, use satellite orbit information, clock information, and inter-code deviation correction models to obtain positioning information and tropospheric characteristic parameter delays .

[0024] (2) The tropospheric zenith dry delay is calculated using the Saastamoinen model using meteorological elements: ; in, is the actual observed value of the tropospheric zenith dry delay, is the geographical latitude of the site, is the air pressure at the measuring station, is the station elevation; (3) Calculate the tropospheric zenith wet delay using the tropospheric zenith total delay and the tropospheric zenith dry delay: ; in, is the actual observed value of the tropospheric zenith wet delay.

[0025] Step 2: Use the historical sounding data of the sounding stations near the station or the reanalysis data of the station location to perform cluster analysis to obtain the statistical characteristics of the regional tropospheric atmospheric parameter profile.

[0026] The clustering algorithm is the Self-Organizing Map (SOM) clustering algorithm. The specific method of using the Self-Organizing Map clustering algorithm to obtain the statistical characteristics of the tropospheric atmospheric parameter profile is as follows: Select the historical sounding data of the sounding station near the station or the reanalysis data of the station location to construct the historical profile data set; pre-process the temperature, relative humidity and refractive index profiles in the data set: first, detect and remove outliers in the original profile data, then calculate the range of meteorological parameter data in each altitude layer, the first-order change rate and second-order change rate of the single profile of the parameter with altitude for each vertical layer data, normalize the data, and unify the parameters and their change rates of different altitude layers; The normalized profile data set is grouped by month, and the data set is subjected to nonlinear dimensionality reduction and topological mapping. The SOM self-organizing mapping network is trained and iteratively adjusted through competitive learning and neighborhood function. For each input data, its Euclidean distance with each neuron node is calculated, and data with similar characteristics are clustered to the same neuron node. After the iteration termination condition is met, the training is completed, and each data is assigned to the cluster represented by the most similar neuron to obtain the final clustering result. For each cluster, its statistical characteristics are calculated. First, the cluster center is extracted, and then the value range of meteorological parameters is hierarchically statistically analyzed along the height dimension. Further, each meteorological element is extracted and its first-order change rate and second-order change rate are calculated to obtain the intra-class statistical distribution characteristics, and finally the statistical characteristics of the tropospheric atmospheric parameter profile are obtained.

[0027] Select grid topology as the topological function to arrange neurons on a two-dimensional grid, and select the distance function as Manhattan distance, which is expressed as: ; in, is the data point, is a randomly initialized neuron weight vector. Find the neuron with the smallest Manhattan distance to the input data point, which is the winning neuron The weight vector is updated according to the winning neuron and the neurons in its neighborhood. The update formula is: ; in, represents the number of iterations, is the learning rate, is the neighborhood function, which represents the neighborhood relationship between the neuron and the winning neuron.

[0028] Step 3: Use the statistical characteristics of the tropospheric atmospheric parameter profile to construct data points on temperature, water vapor pressure and air pressure, and use the ray tracing method or numerical integration method combined with the atmospheric refractive index to forward calculate the theoretical values ​​of the zenith total delay, dry delay and wet delay corresponding to each data point.

[0029] The specific method of constructing data points about temperature, water vapor pressure and air pressure using the statistical characteristics of the tropospheric atmospheric parameter profile is as follows: Assume that the constructed data points are ,in, Indicates temperature, represents the water vapor pressure, Indicates air pressure, Indicates data points, represents the temperature profile, represents the water vapor pressure profile, represents the pressure profile; For temperature profile , is the total number of sections, express The first layer of data Layer; Assume that the minimum value of a single layer of the section is , the maximum value is , based on the profile characteristics, the first-order rate of change of the profile between the extreme values ​​is , whose value range is limited to the minimum value With the maximum value The second-order rate of change is , whose value range is limited to the minimum value With the maximum value between; For the first layer of random numbers , and its generation formula is: ; Here, Representation interval uniformly distributed random variables on , thus ensuring Values ​​in the interval Inside; Subsequent random numbers , for the first-order rate of change limit, by defining the first-order rate of change of each layer , the deformation is: ; in, is the height difference between the upper and lower layers of data; For the second-order rate limit without considering the height variable, it is defined as , the deformation is: ; Take the intersection of the first-order and second-order rate of change constraints for the random numbers of the above two formulas, and finally get the one that satisfies all the constraints. The value range is: ; If the intersection is empty, take the random number generated for the first-order rate of change limit, that is, , Representation interval Uniformly distributed random variables on ; According to the above process, a temperature profile that meets the requirements is generated based on the given single-layer maximum value and the statistical characteristics of the random change rate of the profile. The random number sequence of Sequence and pressure profile sequence, composition as a data point.

[0030] The method of forward modeling using ray tracing is as follows: First, the position of the GNSS antenna at the measuring station is taken as the starting position, and the position vector , ray tracing is performed towards the zenith; the vertical gradient of the atmospheric refractive index is calculated according to the ray tracing step length Interpolate to obtain the atmospheric refractive index corresponding to each calculation step ; From Layer to Layer iterative calculation until the calculation step size The accumulated height reaches the upper limit of the inversion profile; Secondly, the delay of signal propagation is calculated based on the layered atmosphere model. For each segment on the ray path, the refractive index of the segment is used to calculate the delay of signal propagation. The propagation delay of the signal in this section is calculated by the path length, and then the delays of each section are accumulated to obtain the theoretical value of the zenith total delay of the upper limit of the inversion profile; above the upper limit of the profile, the water vapor is zero, and the zenith total delay is represented by the zenith dry delay through the Saastamoinen model; the forward calculation formula for the theoretical value of the total tropospheric propagation delay is as follows: ; in, is the theoretical value of the total delay at the zenith of the troposphere, is the refractive index of each layer of atmosphere, For the ray The path length in the layer, is the zenith dry delay above the upper limit height of the profile; ; in, is the air pressure at the upper limit height of the section, is the latitude at the upper limit of the profile, is the upper limit height of the section; Similarly, the theoretical value of the tropospheric zenith dry delay is calculated as: ; in, is the theoretical value of the tropospheric zenith dry delay, is the dry refractive index of each layer of atmosphere; Theoretical value of tropospheric zenith wet delay The difference between the theoretical value of the total delay at the zenith of the troposphere and the theoretical value of the dry delay at the zenith of the troposphere is obtained: .

[0031] The forward modeling method using numerical integration method is as follows: 1) Assuming that the propagation path of GNSS signals in the troposphere is a straight line, the propagation path of GNSS signals in the troposphere is discretized and divided into several segments. Based on the discretization of the propagation path, the penetration points of each layer are calculated, that is, the points where the GNSS signal propagation path intersects vertically with a specific height layer. The position information of the penetration points is expressed as: ; in, For the The position of the layer's penetration point, is the latitude coordinate of the crossing point, is the longitude coordinate of the crossing point, is the penetration point height; 2) Based on the constructed data points, calculate the atmospheric refractive index profile corresponding to each data point: ; in, Indicates Data point The temperature of the layer, Indicates Data point The water vapor pressure of the layer, Indicates Data point The air pressure of the layer, For the Data point The atmospheric refractive index of the layer is Total Layer, namely: ; 3) The atmospheric refractive index is accumulated along the signal propagation path by the path integration method to obtain the theoretical value of the total delay at the zenith of the troposphere The calculation formula is: ; Theoretical value of tropospheric zenith dry delay The calculation formula is: ; in, is the atmospheric dry refractive index, is the zenith dry delay above the upper limit of the profile, From the Saastamoinen model: ; in, Height of penetration point The air pressure at Height of penetration point The latitude of is the penetration point height; Theoretical value of tropospheric zenith wet delay The theoretical value of the total delay at the zenith of the troposphere Theoretical value of tropospheric zenith dry delay The difference is: .

[0032] In the above scheme, the and The calculation formula is as follows: ; in, is the latitude of the GNSS receiver, is the longitude of the GNSS receiver, For the The satellite azimuth corresponding to the layer-through-layer point, The calculation formula is as follows: ; in, For the The satellite elevation angle corresponding to the layer-through-layer point, is the radius of the Earth, is the GNSS receiver altitude, For the Height of layer penetration point.

[0033] Step 4, select one or more of the total zenith delay, dry delay and wet delay of the troposphere, use the optimization algorithm to invert the tropospheric atmospheric parameter profile, select the data point set with the smallest error between the theoretical value and the actual observation value for iteration, and finally obtain the optimal data point. The temperature, water vapor pressure, air pressure and atmospheric refractive index corresponding to the optimal data point are subjected to data quality control and output as the final inversion result.

[0034] The optimization algorithm includes a single-objective optimization method or a multi-objective optimization method; wherein the single-objective optimization adopts a genetic algorithm or a particle swarm optimization algorithm, and the optimization target is to minimize the root mean square error between the theoretical value and the actual observed value of one of the tropospheric zenith total delay, the tropospheric zenith dry delay and the tropospheric zenith wet delay; The multi-objective optimization adopts non-dominated sorting genetic algorithm or multi-objective particle swarm algorithm, and its optimization goal is to minimize the root mean square error between the theoretical value and the actual observation value of at least two of the tropospheric zenith total delay, tropospheric zenith dry delay and tropospheric zenith wet delay, which can be expressed as follows: ; in, is the multi-objective function vector, is the root mean square error between the actual observed value and the theoretical value of the zenith total delay, is the root mean square error between the actual observed value and the theoretical value of the zenith interference delay, is the root mean square error between the actual observed value and the calculated value of the zenith wet delay, represents the actual observed value of the total delay, Represents the theoretical value of total delay, represents the actual observed value of the dry delay, represents the theoretical value of dry delay, represents the actual observed value of wet delay; Indicates the theoretical wet delay value.

[0035] For stations with co-location observation conditions of ceilometers and microwave radiometers, or conditions for receiving output data from numerical prediction models, the output data of ceilometers, microwave radiometers or numerical prediction models are used as constraints for the inversion algorithm of tropospheric atmospheric parameter profiles; among them, the cloud height data measured by the ceilometer will be used to constrain the relative humidity to be close to saturation near the cloud layer height, with a negative gradient above the cloud layer (decreasing with height) and a positive gradient below the cloud layer (increasing with height); the temperature, air pressure and water vapor pressure profiles output by the microwave radiometer or numerical prediction model are used as the initial values ​​of the data constructed in the subsequent inversion algorithm; If the observation station does not have the conditions for co-location observation of ceilometer and microwave radiometer, and does not have the conditions for receiving output data of numerical prediction model, no constraints will be set.

[0036] When setting constraints, the method of inverting the tropospheric atmospheric parameter profile using the optimization algorithm is as follows: Using Constraint Sets To represent the set of external constraints: Constraints are set using cloud height data measured by ceilometers: profile data is divided into Layers, each layer height is , defines the cloud height measured by the ceilometer Corresponding section layer number , get the layer number closest to the cloud height; according to the existing temperature and air pressure Get saturated water vapor pressure , then the constraint condition for: ; in, Indicates Data point The water vapor pressure of the layer, Indicates Data point The saturated vapor pressure of the layer, is the layer corresponding to the cloud height, Indicates the number of layers corresponding to the cloud height The water vapor pressure, Cloud height corresponding to the number of layers The saturated water vapor pressure; Using the atmospheric parameter profile output by microwave radiometer or numerical prediction model, select the external output profile that can be obtained in actual situation to set the constraint conditions; for the constructed data point profile With the output profile , set the absolute value of the difference between the temperature and water vapor pressure data of each layer of the two to be no greater than a certain minimum value is a constraint condition, then the constraint condition for: ; in, Indicates that the constructed data point is The temperature of the layer, Indicates the external output profile in The temperature of the layer, Indicates that the constructed data point is The water vapor pressure of the layer, Indicates the external output profile in The water vapor pressure of the layer, Indicates that the constructed data point is The air pressure of the layer, Indicates the external output profile in The air pressure of the layer, , and It is The error thresholds of temperature, air pressure and water vapor pressure allowed by the layer; The final objective function is expressed as: ; The point set closest to the actual observation value is selected from the inverted data point set. Through multiple rounds of iterations and gradual optimization, the optimal data points are obtained. The temperature, water vapor pressure, air pressure and atmospheric refractive index corresponding to the optimal data points are subjected to data quality control and then output as the inversion result.

[0037] The data quality control methods are as follows: For multiple solutions, first eliminate the complex solutions; then check the climate limit values, internal consistency of the data and spatial consistency. According to the physical laws and statistical characteristics of the vertical distribution of the troposphere, eliminate the ambiguous solutions that do not conform to the objective physical environment. For example, if the relative humidity exceeds 100%, eliminate the solutions with contradictory relationships between meteorological elements based on the logical relationship between meteorological elements. According to the data of nearby meteorological stations, eliminate the solutions with errors beyond a reasonable range. Finally, the optimal solution of the tropospheric atmospheric parameter profile is obtained.

[0038] The above description of the disclosed embodiments enables one skilled in the art to implement or use the present invention. Various modifications to these embodiments will be apparent to one skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the present invention. Therefore, the present invention will not be limited to the embodiments shown herein, but rather to the widest scope consistent with the principles and novel features disclosed herein.

Claims

1. A method for inverting atmospheric parameter profiles of the marine troposphere based on GNSS, characterized in that: The steps include: Step 1: Use GNSS real-time single-point positioning to estimate the total tropospheric zenith delay of the measuring station, use meteorological element sensors to measure the meteorological elements of the measuring station, calculate the tropospheric zenith dry delay based on the meteorological elements, separate the tropospheric zenith dry delay from the tropospheric zenith total delay, and then obtain the tropospheric zenith wet delay, thereby obtaining the actual observation values ​​of the tropospheric zenith total delay, dry delay and wet delay; Step 2, cluster analysis is performed using historical sounding data from sounding stations near the station or reanalysis data at the station location to obtain statistical characteristics of regional tropospheric atmospheric parameter profiles; Step 3, using the statistical characteristics of the tropospheric atmospheric parameter profile to construct data points on temperature, water vapor pressure and air pressure, and using the ray tracing method or numerical integration method combined with the atmospheric refractive index to forward calculate the theoretical values ​​of the zenith total delay, dry delay and wet delay corresponding to each data point; Step 4, select one or more of the total zenith delay, dry delay and wet delay of the troposphere, use the optimization algorithm to invert the tropospheric atmospheric parameter profile, select the data point set with the smallest error between the theoretical value and the actual observation value for iteration, and finally obtain the optimal data point. The temperature, water vapor pressure, air pressure and atmospheric refractive index corresponding to the optimal data point are subjected to data quality control and output as the final inversion result.

2. The method for inverting the atmospheric parameter profile of the marine troposphere based on GNSS according to claim 1, characterized in that: For stations with co-location observation conditions of ceilometers and microwave radiometers, or conditions for receiving data output from numerical prediction models, the output data of ceilometers, microwave radiometers or numerical prediction models are used as constraints for the inversion algorithm of tropospheric atmospheric parameter profiles; among them, the cloud height data measured by the ceilometer will be used to constrain the relative humidity to be close to saturation near the cloud layer height, with a negative gradient above the cloud layer and a positive gradient below the cloud layer; the temperature, air pressure and water vapor pressure profiles output by the microwave radiometer or numerical prediction model are used as the initial values ​​of the data constructed in the subsequent inversion algorithm; If the observation station does not have the conditions for co-location observation of ceilometer and microwave radiometer, and does not have the conditions for receiving output data of numerical prediction model, no constraints will be set.

3. The method for inverting the atmospheric parameter profile of the marine troposphere based on GNSS according to claim 1, characterized in that: In step 1, the method for obtaining the actual observed values ​​of the total tropospheric zenith delay, dry delay, and wet delay is as follows: (1) Using GNSS receivers to obtain dual-frequency raw observations, combined with Beidou satellite-based broadcast PPP-B2b corrections, the ionospheric delay error is eliminated through ionospheric-free combination, and the actual observation value of the total tropospheric zenith delay is extracted. ; (2) Calculate the tropospheric zenith dry delay using meteorological elements: ； in, is the actual observed value of the tropospheric zenith dry delay, is the geographical latitude of the site, is the air pressure at the measuring station, is the station elevation; (3) Calculate the tropospheric zenith wet delay using the tropospheric zenith total delay and the tropospheric zenith dry delay: ； in, is the actual observed value of the tropospheric zenith wet delay.

4. The method for inverting the atmospheric parameter profile of the marine troposphere based on GNSS according to claim 1, characterized in that: The clustering algorithm is a self-organizing map clustering algorithm. The specific method of using the self-organizing map clustering algorithm to obtain the statistical characteristics of the tropospheric atmospheric parameter profile is as follows: Select the historical sounding data of the sounding station near the station or the reanalysis data of the station location to construct the historical profile data set; pre-process the temperature, relative humidity and refractive index profiles in the data set: first, detect and remove outliers in the original profile data, then calculate the range of meteorological parameter data in each altitude layer, the first-order change rate and second-order change rate of the single profile of the parameter with altitude for each vertical layer data, normalize the data, and unify the parameters and their change rates of different altitude layers; The normalized profile data sets were grouped by month, nonlinear dimension reduction and topological mapping were performed on the data sets, and the SOM self-organizing map network was trained through competitive learning and iterative adjustment of neighborhood functions. For each input data, the Euclidean distance between it and each neuron node is calculated, and the data with similar characteristics are clustered to the same neuron node. After the iteration termination condition is met, the training is completed, and each data is assigned to the cluster represented by the neuron that is most similar to it to obtain the final clustering result; for each cluster, its statistical characteristics are calculated. First, the cluster center is extracted, and then the value range of the meteorological parameters is hierarchically statistically analyzed along the height dimension. Further, each meteorological element is extracted to calculate its first-order change rate and second-order change rate, and the statistical distribution characteristics within the class are obtained. Finally, the statistical characteristics of the tropospheric atmospheric parameter profile are obtained.

5. The method for inverting the atmospheric parameter profile of the marine troposphere based on GNSS according to claim 1, characterized in that: In step 3, the specific method of constructing data points about temperature, water vapor pressure and air pressure by using the statistical characteristics of the tropospheric atmospheric parameter profile is as follows: Assume that the constructed data points are ,in, Indicates temperature, represents the water vapor pressure, Indicates air pressure, Indicates data points, represents the temperature profile, represents the water vapor pressure profile, represents the pressure profile; For temperature profile , is the total number of sections, express The first layer of data Layer; Assume that the minimum value of a single layer of the section is , the maximum value is , based on the profile characteristics, the first-order rate of change of the profile between the extreme values ​​is , whose value range is limited to the minimum value With the maximum value The second-order rate of change is , whose value range is limited to the minimum value With the maximum value between; For the first layer of random numbers , and its generation formula is: ； Here, Representation interval uniformly distributed random variables on , thus ensuring Values ​​in the interval Inside; Subsequent random numbers , for the first-order rate of change limit, by defining the first-order rate of change of each layer , the deformation is: ； in, is the height difference between the upper and lower layers of data; For the second-order rate limit without considering the height variable, it is defined as , the deformation is: ； Take the intersection of the first-order and second-order rate of change constraints for the random numbers of the above two formulas, and finally get the one that satisfies all the constraints. The value range is: ； If the intersection is empty, take the random number generated for the first-order rate of change limit, that is, , Representation interval Uniformly distributed random variables on ; According to the above process, a temperature profile that meets the requirements is generated based on the given single-layer maximum value and the statistical characteristics of the random change rate of the profile. The random number sequence of Sequence and pressure profile sequence, composition as a data point.

6. The method for inverting the atmospheric parameter profile of the marine troposphere based on GNSS according to claim 1, characterized in that: In step 3, the method of forward modeling using ray tracing is as follows: First, the position of the GNSS antenna at the measuring station is taken as the starting position, and the position vector , ray tracing is performed towards the zenith; Calculate the vertical gradient of atmospheric refractivity by ray tracing Interpolate to obtain the atmospheric refractive index corresponding to each calculation step ; From Layer to Layer iterative calculation until the calculation step size The accumulated height reaches the upper limit of the inversion profile; Secondly, the delay of signal propagation is calculated based on the layered atmosphere model. For each segment on the ray path, the refractive index of the segment is used to calculate the delay of signal propagation. The propagation delay of the signal in this section is calculated by the path length, and then the delays of each section are accumulated to obtain the theoretical value of the zenith total delay of the upper limit of the inversion profile; above the upper limit of the profile, the water vapor is zero, and the zenith total delay is represented by the zenith dry delay through the Saastamoinen model; the forward calculation formula for the theoretical value of the total tropospheric propagation delay is as follows: ； in, is the theoretical value of the total delay at the zenith of the troposphere, is the refractive index of each layer of atmosphere, For the ray The path length in the layer, is the zenith dry delay above the upper limit height of the profile; ； in, is the air pressure at the upper limit height of the section, is the latitude at the upper limit of the profile, is the upper limit height of the section; Similarly, the theoretical value of the tropospheric zenith dry delay is calculated as: ； in, is the theoretical value of the tropospheric zenith dry delay, is the dry refractive index of each layer of atmosphere; Theoretical value of tropospheric zenith wet delay The difference between the theoretical value of the total delay at the zenith of the troposphere and the theoretical value of the dry delay at the zenith of the troposphere is obtained: 。 7. The method for inverting the atmospheric parameter profile of the marine troposphere based on GNSS according to claim 1, characterized in that: In step 3, the forward modeling method using the numerical integration method is as follows: 1) Assuming that the propagation path of GNSS signals in the troposphere is a straight line, the propagation path of GNSS signals in the troposphere is discretized and divided into several segments. Based on the discretization of the propagation path, the penetration points of each layer are calculated, that is, the points where the GNSS signal propagation path intersects vertically with a specific height layer. The position information of the penetration points is expressed as: ； in, For the The position of the layer's penetration point, is the latitude coordinate of the crossing point, is the longitude coordinate of the crossing point, is the penetration point height; 2) Based on the constructed data points, calculate the atmospheric refractive index profile corresponding to each data point: ； in, Indicates Data point The temperature of the layer, Indicates Data point The water vapor pressure of the layer, Indicates Data point The air pressure of the layer, For the Data point The atmospheric refractive index of the layer is Total Layer, namely: ； 3) The atmospheric refractive index is accumulated along the signal propagation path by the path integration method to obtain the theoretical value of the total delay at the zenith of the troposphere The calculation formula is: ； Theoretical value of tropospheric zenith dry delay The calculation formula is: ； in, is the atmospheric dry refractive index, is the zenith dry delay above the upper limit of the profile, From the Saastamoinen model: ； in, Height of penetration point The air pressure at Height of penetration point The latitude of is the penetration point height; Theoretical value of tropospheric zenith wet delay The theoretical value of the total delay at the zenith of the troposphere Theoretical value of tropospheric zenith dry delay The difference is obtained: 。 8. The method for inverting the atmospheric parameter profile of the marine troposphere based on GNSS according to claim 7, characterized in that: Each layer point and The calculation formula is as follows: ； in, is the latitude of the GNSS receiver, is the longitude of the GNSS receiver, For the The satellite azimuth corresponding to the layer-through-layer point, The calculation formula is as follows: ； in, For the The satellite elevation angle corresponding to the layer-through-layer point, is the radius of the Earth, is the GNSS receiver altitude, For the Height of layer penetration point.

9. The method for inverting the atmospheric parameter profile of the marine troposphere based on GNSS according to claim 1, characterized in that: In step 4, the optimization algorithm includes a single-objective optimization method or a multi-objective optimization method; wherein the single-objective optimization adopts a genetic algorithm or a particle swarm optimization algorithm, and the optimization goal is to minimize the root mean square error between the theoretical value and the actual observed value of one of the tropospheric zenith total delay, the tropospheric zenith dry delay, and the tropospheric zenith wet delay; The multi-objective optimization adopts non-dominated sorting genetic algorithm or multi-objective particle swarm algorithm, and its optimization goal is to minimize the root mean square error between the theoretical value and the actual observation value of at least two of the tropospheric zenith total delay, tropospheric zenith dry delay and tropospheric zenith wet delay, which can be expressed as follows: ； in, is the multi-objective function vector, is the root mean square error between the actual observed value and the theoretical value of the zenith total delay, is the root mean square error between the actual observed value and the theoretical value of the zenith interference delay, is the root mean square error between the actual observed value and the calculated value of the zenith wet delay, represents the actual observed value of the total delay, Represents the theoretical value of total delay, represents the actual observed value of the dry delay, represents the theoretical value of dry delay, represents the actual observed value of wet delay; Indicates the theoretical wet delay value.

10. The method for inverting the atmospheric parameter profile of the marine troposphere based on GNSS according to claim 2, characterized in that: In step 4, when setting constraints, the method of inverting the tropospheric atmospheric parameter profile using the optimization algorithm is as follows: Using Constraint Sets To represent the set of external constraints: Constraints are set using cloud height data measured by ceilometers: profile data is divided into Layers, each layer height is , defines the cloud height measured by the ceilometer Corresponding section layer number , get the layer number closest to the cloud height; according to the existing temperature and air pressure Get saturated water vapor pressure , then the constraint condition for: ； in, Indicates Data point The water vapor pressure of the layer, Indicates Data point The saturated vapor pressure of the layer, is the layer corresponding to the cloud height, Indicates the number of layers corresponding to the cloud height The water vapor pressure, Indicates the number of layers corresponding to the cloud height saturated water vapor pressure; Using the atmospheric parameter profile output by microwave radiometer or numerical prediction model, select the external output profile that can be obtained in actual situation to set the constraint conditions; for the constructed data point profile With the output profile , set the absolute value of the difference between the temperature and water vapor pressure data of each layer of the two to be no greater than a certain minimum value is a constraint condition, then the constraint condition for: ； in, Indicates that the constructed data point is The temperature of the layer, Indicates the external output profile in The temperature of the layer, Indicates that the constructed data point is The water vapor pressure of the layer, Indicates the external output profile in The water vapor pressure of the layer, Indicates that the constructed data point is The air pressure of the layer, Indicates the external output profile in The air pressure of the layer, , and It is The error thresholds of temperature, air pressure and water vapor pressure allowed by the layer; The final objective function is expressed as: ； The point set closest to the actual observation value is selected from the inverted data point set. Through multiple rounds of iterations and gradual optimization, the optimal data points are obtained. The temperature, water vapor pressure, air pressure and atmospheric refractive index corresponding to the optimal data points are subjected to data quality control and then output as the inversion result.

Images