A method for retrieving the tropospheric atmospheric parameter profile over the sea based on GNSS

Abstract

The present invention relates to the field of marine meteorological monitoring, and discloses a method for retrieving the tropospheric atmospheric parameter profile based on GNSS, which includes the following steps: obtaining the actual observed values of the total zenith delay, dry delay, and wet delay of the troposphere; obtaining the statistical characteristics of the tropospheric atmospheric parameter profile through cluster analysis; constructing data points about temperature, water vapor pressure, and air pressure by using the statistical characteristics of the tropospheric atmospheric parameter profile, and forward calculating the theoretical values of the total zenith delay, dry delay, and wet delay corresponding to each data point; selecting one or more of the total zenith delay, dry delay, and wet delay of the troposphere, and using an optimization algorithm to retrieve the tropospheric atmospheric parameter profile, iteratively obtaining the optimal data points, and outputting the temperature, water vapor pressure, air pressure, and atmospheric refractive index corresponding to the optimal data points after data quality control. The method disclosed by the present invention has a relatively low equipment installation and maintenance cost, is suitable for marine observation platforms, has reliable inversion results, and can provide support for marine tropospheric atmospheric research.

Description

Technical Field

[0001] The present invention relates to the field of marine meteorological monitoring, and particularly to a method for retrieving the tropospheric atmospheric parameter profile based on GNSS in the sea area. Background Art

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

[0003] The methods for obtaining the tropospheric atmospheric parameter profile are divided into two types: contact type and non-contact type. The contact type usually uses the method of releasing radiosondes and sounding rockets. The non-contact methods include microwave radiometers, lidars, etc. The above detection methods have certain deficiencies when applied to marine observation platforms. The time resolution of releasing radiosondes and sounding rockets is limited, and manual release is required, which has low operability on marine platforms and cannot achieve continuous marine observation. Microwave radiometers can achieve a high 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, it is difficult to be applied to sea-based observation platforms. Detection by optical methods such as lidars is easily affected by the high-humidity and high-salt environment at sea and is not corrosion-resistant. Although the technology of autonomous sounding observation using unmanned aerial vehicles and unmanned ships has developed rapidly in recent years, in actual use, the unmanned sounding system in the marine environment needs to solve the problem of replenishing sounding sensors. Summary of the Invention

[0004] To solve the above technical problems, the present invention provides a method for retrieving the tropospheric atmospheric parameter profile based on GNSS in the sea area, so as to achieve the purpose of low equipment installation and maintenance costs, being suitable for marine observation platforms, having reliable inversion results, and providing support for marine tropospheric atmosphere research.

[0005] To achieve the above object, the technical solution of the present invention is as follows:

[0006] A method for retrieving the tropospheric atmospheric parameter profile based on GNSS in the sea area includes the following steps:

[0007] Step 1: Use GNSS real-time single-point positioning to estimate the total zenith tropospheric delay of the station position, measure the meteorological elements at the station position using a meteorological element sensor, calculate the dry zenith tropospheric delay based on the meteorological elements, separate the dry zenith tropospheric delay from the total zenith tropospheric delay to obtain the wet zenith tropospheric delay, so as to obtain the actual observed values of the total zenith tropospheric delay, dry delay, and wet delay;

[0008] Step 2: Use the historical sounding data of the sounding station near the station or the reanalysis data of the station position for cluster analysis to obtain the statistical characteristics of the regional tropospheric atmospheric parameter profile;

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

[0010] Step 4: Select one or more of the total zenith tropospheric delay, dry delay, and wet delay, use an optimization algorithm to invert the tropospheric atmospheric parameter profile, select the data point set with a smaller error between the theoretical value and the actual observed value for iteration, and finally obtain the optimal data points. After performing data quality control on the temperature, water vapor pressure, air pressure, and atmospheric refractive index corresponding to the optimal data points, output them as the final inversion result.

[0011] In the above solution, for a station with the co-location observation conditions of a ceilometer and a microwave radiometer, or the data reception conditions of a numerical weather prediction model output, use the ceilometer, microwave radiometer, or numerical weather prediction model output data as the constraint conditions for 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 near saturation at the cloud 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 weather prediction model are used as the initial values of the data constructed in the subsequent inversion algorithm;

[0012] If the station does not have the co-location observation conditions of a ceilometer and a microwave radiometer, as well as the data reception conditions of a numerical weather prediction model output, no constraint conditions will be set.

[0013] In the above solution, in Step 1, the method for obtaining the actual observed values of the total zenith tropospheric delay, dry delay, and wet delay is as follows:

[0014] (1) Use a GNSS receiver to obtain dual-frequency raw observations, combine the Beidou satellite-based broadcast PPP-B2b corrections, eliminate the ionospheric delay error through the ionosphere-free combination, and extract the actual observed value of the total zenith tropospheric delay ;

[0015] (2) Calculate the dry zenith tropospheric delay using meteorological elements:

[0016] ;

[0017] Among them, is the actual observed value of the tropospheric zenith dry delay, is the geographical latitude of the station, is the air pressure at the station location, is the elevation of the station;

[0018] (3) Calculate the tropospheric zenith wet delay by using the total tropospheric zenith delay and the tropospheric zenith dry delay:

[0019] ;

[0020] Among them, is the actual observed value of the tropospheric zenith wet delay.

[0021] In the above scheme, the clustering algorithm is the self-organizing mapping clustering algorithm. The specific method for obtaining the statistical characteristics of the tropospheric atmospheric parameter profile by using the self-organizing mapping clustering algorithm is as follows:

[0022] Select the historical sounding data of the sounding station near the station or the reanalysis data of the station location to construct a historical profile data set; preprocess the temperature, relative humidity, and refractive index profiles in the data set: first, detect and remove outliers from the original profile data, and then, for the vertical stratification data of each profile, calculate the value range of the meteorological parameter data, the first-order change rate and the second-order change rate of the parameter with height layer by layer, normalize the data, and unify the parameters and their change rates in different height layers;

[0023] Group the normalized profile data set by month, perform nonlinear dimensionality reduction and topological mapping on the data set, train the SOM self-organizing mapping network, and iteratively adjust through competitive learning and neighborhood function; for each input data, calculate its Euclidean distance from each neuron node, cluster the data with similar characteristics into the same neuron node, and after meeting the iterative termination condition, the training is completed, and each data is assigned to the cluster represented by the neuron most similar to it to obtain the final clustering result; calculate the statistical characteristics for each cluster, first extract the cluster center, secondly statistically analyze the value range of the meteorological parameters layer by layer along the height dimension, and further extract the first-order change rate and the second-order change rate of each meteorological element to obtain the intra-class statistical distribution characteristics, and finally obtain the statistical characteristics of the tropospheric atmospheric parameter profile.

[0024] In the above scheme, in step 3, the specific method for 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:

[0025] Assume that the constructed data point is , where, represents temperature, represents water vapor pressure, represents air pressure, represents the th data point, represents the temperature profile, represents the water vapor pressure profile, represents the air pressure profile;

[0026] For the temperature profile , is the total number of layers in the profile, represents the th layer in the layer data; assuming the minimum value of a single layer in the profile is and the maximum value is , the first-order change rate of the profile limited by the minimum and maximum values of the profile is , and its value range is defined between the minimum value and the maximum value , and the second-order change rate is , and its value range is defined between the minimum value and the maximum value ;

[0027] For the random number of the first layer, its generation formula is:

[0028] ;

[0029] Here, represents a uniformly distributed random variable on the interval to ensure that takes values within the interval ;

[0030] For subsequent random numbers , for the first-order change rate limit, by defining the first-order change rate of each layer, it is transformed to get:

[0031] ;

[0032] where, is the height difference between the data of the upper and lower layers;

[0033] For the second-order change rate limit without considering the height variable, it is defined as , and it is transformed to get:

[0034] ;

[0035] Take the intersection of the first-order and second-order change rate limit conditions of the random numbers in the above two formulas, and finally obtain the value range is:

[0036] ;

[0037] If the intersection is empty, then take the random numbers generated for the first-order change rate limit, that is , represents a uniformly distributed random variable on the interval ;

[0038] According to the above process, based on the statistical characteristics of the single-layer maximum and minimum values of the given profile and the random change rate, generate a random number sequence that meets the requirements for the temperature profile , and similarly generate a water vapor pressure profile sequence and a barometric pressure profile sequence, which constitutes as data points.

[0039] In the above solution, in step 3, the method of forward modeling using the ray tracing method is as follows:

[0040] First, take the position of the GNSS antenna at the measuring station as the starting position, and the position vector , and perform ray tracing in the zenith direction; interpolate the vertical gradient of the atmospheric refractive index according to the ray tracing calculation step size to obtain the atmospheric refractive index corresponding to each calculation step size; iterate from the th layer to the th layer until the cumulative height of the calculation step size reaches the upper limit of the inversion profile;

[0041] Secondly, calculate the delay of signal propagation based on the stratified atmosphere model. For each segment on the ray path, calculate the propagation delay of the signal in this segment according to the refractive index of this segment and the path length, and then accumulate the delays of each segment to obtain the theoretical value of the total zenith delay at the upper limit of the inversion profile; above the upper limit height of the profile, the water vapor is zero, and the total zenith delay is represented by the Saastamoinen model by the zenith dry delay; then the formula for the theoretical value of the total tropospheric propagation delay in forward modeling is as follows:

[0042] ;

[0043] Among them, is the theoretical value of the total tropospheric zenith delay, is the atmospheric refractive index of each layer, is the path length of the ray in the th layer, is the zenith dry delay above the upper limit height of the section;

[0044] ;

[0045] wherein, is the air pressure at the upper limit height of the section, is the latitude at the upper limit height of the section, is the upper limit height of the section;

[0046] Similarly, calculate the theoretical value of the tropospheric zenith dry delay:

[0047] ;

[0048] wherein, is the theoretical value of the tropospheric zenith dry delay, is the dry refractivity of each layer of the atmosphere;

[0049] The theoretical value of the tropospheric zenith wet delay is obtained by subtracting the theoretical value of the tropospheric zenith dry delay from the theoretical value of the tropospheric zenith total delay:

[0050] .

[0051] In the above solution, in step 3, the method of forward modeling using the numerical integration method is as follows:

[0052] 1) Assume that the propagation path of the GNSS signal in the troposphere is a straight line, discretize the propagation path of the GNSS signal in the troposphere, divide the propagation path into several segments, and on the basis of the discretization of the propagation path, calculate the penetration points of each layer, that is, the points where the GNSS signal propagation path intersects perpendicularly with a specific height level. The position information of the penetration points is expressed as:

[0053] ;

[0054] wherein, is the position of the penetration point of the th layer, is the latitude coordinate of the penetration point, is the longitude coordinate of the penetration point, is the height of the penetration point;

[0055] 2) According to the constructed data points, calculate the atmospheric refractive index profile corresponding to each data point:

[0056] ;

[0057] wherein, represents the temperature of the th data point in the th layer, represents the The water vapor pressure at the layer, represents the atmospheric pressure at the layer, is the atmospheric refractive index at the layer, then there are a total of layers, namely:

[0058] ;

[0059] 3) By using the path integral method, the atmospheric refractive index is accumulated along the signal propagation path to obtain the theoretical value of the tropospheric zenith total delay The calculation formula is:

[0060] ;

[0061] The theoretical value of the tropospheric zenith dry delay The calculation formula is:

[0062] ;

[0063] Among them, is the atmospheric dry refractive index, is the zenith dry delay above the upper limit height of the profile, obtained from the Saastamoinen model:

[0064] ;

[0065] Among them, is the height of the layer crossing point where the atmospheric pressure is is the latitude at the height of the layer crossing point where the height is ;

[0066] The theoretical value of the tropospheric zenith wet delay is obtained by subtracting the theoretical value of the tropospheric zenith dry delay from the theoretical value of the tropospheric zenith total delay :

[0067] .

[0068] In the above scheme, the and of each layer crossing point are calculated as follows:

[0069] ;

[0070] Among them, is the latitude of the GNSS receiver, is the longitude of the GNSS receiver, is the satellite azimuth angle corresponding to the penetration point of the layer, and its calculation formula is as follows:

[0071] ;

[0072] wherein, is the satellite elevation angle corresponding to the penetration point of the layer, is the radius of the earth, is the height of the GNSS receiver, is the height of the penetration point of the

[0073] In the above solution, in step 4, the optimization algorithm includes a single-objective optimization or a multi-objective optimization method; wherein, the single-objective optimization adopts a genetic algorithm or a particle swarm optimization algorithm, and the optimization objective 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;

[0074] The multi-objective optimization adopts a non-dominated sorting genetic algorithm or a multi-objective particle swarm algorithm, and its optimization objective is to minimize the root mean square error value between the theoretical value and the actual observed value of at least two of the tropospheric zenith total delay, the tropospheric zenith dry delay, and the tropospheric zenith wet delay, which is expressed by the following formula:

[0075] ;

[0076] wherein, is the multi-objective function vector, is the root mean square error value between the actual observed value and the theoretical value of the zenith total delay, is the root mean square error value between the actual observed value and the theoretical value of the zenith dry delay, is the root mean square error value 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 the total delay, represents the actual observed value of the dry delay, represents the theoretical value of the dry delay, represents the actual observed value of the wet delay; represents the theoretical value of the wet delay.

[0077] In a further technical solution, in step 4, when setting the constraint conditions, the method for inversely calculating the tropospheric atmospheric parameter profile using the optimization algorithm is as follows:

[0078] Use the constraint set is used to represent the set of external constraint conditions:

[0079] Set the constraint conditions using the cloud height data measured by the ceilometer: The profile data is divided into layers, and the height of each layer is . Define the cloud height measured by the ceilometer and the corresponding profile layer number to obtain the layer number closest to the cloud height; According to the existing temperature and air pressure to obtain the saturated water vapor pressure . Then the constraint condition is:

[0080] ;

[0081] Among them, represents the water vapor pressure of the st data point at the th layer, represents the saturated water vapor pressure of the st data point at the th layer, is the layer corresponding to the cloud height, represents the water vapor pressure of the layer corresponding to the cloud height , the saturated water vapor pressure of the layer corresponding to the cloud height ;

[0082] Use the atmospheric parameter profiles output by the microwave radiometer or the numerical weather prediction model, and select the external output profile that can be obtained in the actual situation among the two to set the constraint conditions; For the constructed data point profile and 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 small value as the constraint condition. Then the constraint condition is:

[0083] ;

[0084] Among them, represents the temperature of the constructed data point at the th layer, represents the temperature of the external output profile at the th layer, represents the water vapor pressure of the constructed data point at the th layer, represents the water vapor pressure of the external output profile at the th layer, represents the air pressure of the constructed data point at the th layer, represents the external output profile at the The air pressure of the layer, 、 and are the error thresholds of the allowable temperature, air pressure and water vapor pressure of the th layer;

[0085] Then the final objective function is expressed as:

[0086] ;

[0087] Select the point set closest to the actual observation value from the data point set obtained by inversion. Through multiple rounds of iteration and gradual optimization, the optimal data points are obtained. After data quality control of the temperature, water vapor pressure, air pressure and atmospheric refractivity corresponding to the optimal data points, they are output as the inversion result.

[0088] Through the above technical solutions, a method for inverting the tropospheric atmospheric parameter profile based on GNSS provided by the present invention has the following beneficial effects:

[0089] 1. Through the collaborative observation of GNSS signals and meteorological sensors, combined with the Beidou satellite-based broadcast PPP-B2b correction and dual-frequency ionospheric error elimination technology, and using the Saastamoinen model, the present invention obtains high-precision actual observation values of the total zenith delay, dry delay and wet delay of the troposphere;

[0090] 2. The present invention adopts the self-organizing mapping (SOM) clustering algorithm to construct a statistical characteristic model of the tropospheric atmospheric parameter profile, and constructs an initial value for subsequent inversion, thereby improving the inversion accuracy;

[0091] 3. The present invention uses the 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 zenith delay, dry delay and wet delay of the troposphere, providing accurate theoretical data for setting the objective function for subsequent inversion;

[0092] 4. The present invention performs inversion through an optimization algorithm combined with dynamic external constraint conditions, considers using a multi-objective function as the objective, and sets external constraint conditions for the inversion algorithm according to specific situations, further improving the accuracy of the inversion result, thereby obtaining the optimal solution. In summary, the method for inverting the tropospheric atmospheric parameter profile based on GNSS signals provided by the present invention can realize the inversion of the tropospheric atmospheric parameter profile without relying on historical sounding data, has the unique advantage of being applicable to marine observation platforms lacking sounding data, has a relatively low equipment installation and maintenance cost, and can provide support for marine tropospheric atmospheric research. BRIEF DESCRIPTION OF THE DRAWINGS

[0093] In order to more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the following will briefly introduce the drawings required for the description of the embodiments or the prior art.

[0094] Figure 1 Schematic diagram of the process of a method for retrieving tropospheric atmospheric parameter profiles based on GNSS disclosed in an embodiment of the present invention;

[0095] Figure 2 Schematic diagram of the installation of the station equipment involved in the present invention. Detailed implementation manners

[0096] Next, the technical solutions in the embodiments of the present invention will be clearly and completely described in conjunction with the accompanying drawings in the embodiments of the present invention.

[0097] The present invention provides a method for retrieving tropospheric atmospheric parameter profiles based on GNSS, as Figure 1 shown, and includes the following steps:

[0098] Step 1, use GNSS real-time single-point positioning to estimate the tropospheric zenith total delay of the station position, measure the meteorological elements at the station position using a meteorological element sensor, calculate the tropospheric zenith dry delay according to the meteorological elements, separate the tropospheric zenith dry delay from the tropospheric zenith total delay to obtain the tropospheric zenith wet delay, so as to obtain the actual observed values of the tropospheric zenith total delay, dry delay and wet delay.

[0099] The installation of the station equipment involved in the present invention is as Figure 2 shown, and the main equipment involved includes: external equipment such as a meteorological element sensor, a GNSS antenna, a GNSS receiver, a data processing unit, and a ceilometer. Among them, the meteorological element sensor and the GNSS antenna 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 situation of the installation platform. A method for retrieving tropospheric atmospheric parameter profiles based on GNSS provided by the present invention, the input data of which is collected from a meteorological element sensor and a GNSS receiver, and the algorithm software of this method runs on the data processing unit. For a station that already has a meteorological element sensor and a GNSS receiver, only a data processing unit and the corresponding algorithm software need to be added to conveniently implement the function of retrieving tropospheric atmospheric parameter profiles, and the equipment installation and maintenance costs are relatively low.

[0100] The method for obtaining the actual observed values of the tropospheric zenith total delay, dry delay and wet delay is as follows:

[0101] (1) Use a GNSS receiver to obtain dual-frequency raw observations, combine the Beidou satellite-based broadcast PPP-B2b corrections, eliminate the ionospheric delay error through the ionosphere-free combination, and extract the actual observed value of the tropospheric zenith total delay ;

[0102] The air pressure at the GNSS antenna position 、temperature and relative humidity and other meteorological elements are measured by meteorological element sensors; the GNSS raw observables pseudorange , carrier phase , Doppler frequency shift and carrier-to-noise ratio are obtained by the GNSS receiver, and the corrections are obtained by the GNSS receiver or an external computing unit.

[0103] The signal transmitted by the GNSS antenna is processed by the signal processing board, converted into a digital signal, and the satellite signal is searched through correlation operations in the digital domain to complete signal acquisition. The navigation data in the signal is demodulated, and the binary data sequence containing parameters such as satellite orbit, clock, and ionospheric model is extracted by operating in reverse according to the modulation method and transmitted to the data processing module.

[0104] Real-time precise point positioning is performed using the high-precision GNSS data processing software RTKLIB. Receive and parse the observation data and navigation message from the GNSS receiver, and perform quality control and preprocessing on the pseudorange and carrier phase information in the observation data. Through precise point positioning, use satellite orbit information, clock error information, and inter-code bias correction model to obtain positioning information and tropospheric characteristic parameter delay .

[0105] (2) Calculate the tropospheric zenith dry delay using the Saastamoinen model with meteorological elements:

[0106] ;

[0107] where is the actual observed value of the tropospheric zenith dry delay, is the geographical latitude of the station, is the air pressure at the station location, is the elevation of the station;

[0108] (3) Calculate the tropospheric zenith wet delay using the tropospheric zenith total delay and the tropospheric zenith dry delay:

[0109] ;

[0110] where is the actual observed value of the tropospheric zenith wet delay.

[0111] Step 2: Perform cluster analysis using the historical sounding data of the sounding station near the station or the reanalysis data of the station location to obtain the statistical characteristics of the regional tropospheric atmospheric parameter profile.

[0112] The clustering algorithm is the Self-Organizing Map (SOM) clustering algorithm. The specific method for obtaining the statistical characteristics of the tropospheric atmospheric parameter profiles using the self-organizing map clustering algorithm is as follows:

[0113] Select the historical sounding data of the sounding station near the measurement station or the reanalysis data of the measurement station location to construct a historical profile dataset; preprocess the temperature, relative humidity, and refractive index profiles in the dataset: first, detect and eliminate outliers from the original profile data, and then calculate the value range of the meteorological parameter data, the first-order change rate and the second-order change rate of the parameter with height layer by layer for the vertical stratification data of each profile, normalize the data, and unify the parameters and their change rates at different height layers;

[0114] Group the normalized profile dataset by month, perform non-linear dimensionality reduction and topological mapping on the dataset, train the SOM self-organizing map network, and iteratively adjust through competitive learning and neighborhood functions; for each input data, calculate its Euclidean distance from each neuron node, cluster the data with similar characteristics into the same neuron node, and after meeting the iteration termination condition, the training is completed, and each data is assigned to the cluster represented by the neuron most similar to it to obtain the final clustering result; calculate the statistical characteristics for each cluster, first extract the cluster center, then statistically calculate the value range of the meteorological parameters layer by layer along the height dimension, further extract each meteorological element and calculate its first-order change rate and second-order change rate to obtain the intra-class statistical distribution characteristics, and finally obtain the statistical characteristics of the tropospheric atmospheric parameter profiles.

[0115] Select the grid topology as the topological function to arrange the neurons on a two-dimensional grid, and select the distance function as the Manhattan distance, and its expression is:

[0116] ;

[0117] where is the data point, is the randomly initialized neuron weight vector. Find the neuron with the minimum Manhattan distance from the input data point, which is the winning neuron . Update the weight vector according to the winning neuron and the neurons within its neighborhood, and the update formula is:

[0118] ;

[0119] where represents the number of iterations, is the learning rate, is the neighborhood function, indicating the neighborhood relationship between the neuron and the winning neuron.

[0120] Step 3: Construct data points regarding temperature, water vapor pressure, and air pressure by using the statistical characteristics of the tropospheric atmospheric parameter profile. Forward calculate the theoretical values of the zenith total delay, dry delay, and wet delay corresponding to each data point by using the ray tracing method or numerical integration method in combination with the atmospheric refractive index.

[0121] The specific method for constructing data points regarding temperature, water vapor pressure, and air pressure by using the statistical characteristics of the tropospheric atmospheric parameter profile is as follows:

[0122] Assume that the constructed data point is , where represents temperature, represents water vapor pressure, represents air pressure, represents the th data point, represents the temperature profile, represents the water vapor pressure profile, represents the air pressure profile;

[0123] For the temperature profile , is the total number of layers in the profile, represents the th layer in the layer data; assume that the minimum value of a single layer in the profile is , the maximum value is , the first-order change rate of the profile limited between the minimum and maximum values is and its value range is defined between the minimum value and the maximum value , the second-order change rate is and its value range is defined between the minimum value and the maximum value;

[0124] For the first-layer random number , its generation formula is:

[0125] ;

[0126] Here, represents a uniformly distributed random variable on the interval to ensure that takes values within the interval ;

[0127] For subsequent random numbers , for the first-order change rate limit, by defining the first-order change rate of each layer, it is transformed to get:

[0128] ;

[0129] Among them, is the height difference between the upper and lower layers of data;

[0130] For the second-order change rate limit without considering the height variable, it is defined as , and by transformation, we get:

[0131] ;

[0132] Take the intersection of the first-order and second-order change rate limit conditions for the random numbers of the above two formulas, and finally obtain the value range that satisfies all the limit conditions is:

[0133] ;

[0134] If the intersection is empty, then take the random number generated for the first-order change rate limit, that is , represents a uniformly distributed random variable on the interval ;

[0135] According to the above process, based on the statistical characteristics of the given single-layer maximum value and random change rate of the profile, generate a random number sequence that meets the requirements for the temperature profile , and similarly generate a sequence for the water vapor pressure profile and a sequence for the air pressure profile , which constitute as data points.

[0136] The method of forward modeling using the ray tracing method is as follows:

[0137] First, take the position of the GNSS antenna at the station as the starting position, with the position vector , and perform ray tracing in the zenith direction; interpolate the vertical gradient of the atmospheric refractive index according to the ray tracing calculation step size to obtain the atmospheric refractive index corresponding to each calculation step size; iterate from the th layer to the th layer until the cumulative height of the calculation step size reaches the upper limit of the inversion profile;

[0138] Secondly, calculate the delay of signal propagation based on the stratified atmospheric model. For each segment on the ray path, calculate the propagation delay of the signal in this segment according to the refractive index of this segment and the path length, and then accumulate the delays of each segment to obtain the theoretical value of the total zenith delay at the upper limit of the inversion profile; above the upper limit height of the profile, the water vapor is zero, and the total zenith delay is represented by the Saastamoinen model through the zenith dry delay; then the theoretical value formula for calculating the total tropospheric propagation delay in forward modeling is as follows:

[0139] ;

[0140] Among them, is the theoretical value of the tropospheric zenith total delay, is the atmospheric refractive index of each layer, is the path length of the ray in the th layer, is the zenith dry delay above the upper limit height of the profile;

[0141] ;

[0142] Among them, is the air pressure at the upper limit height of the profile, is the latitude at the upper limit height of the profile, is the upper limit height of the profile;

[0143] Similarly, calculate the theoretical value of the tropospheric zenith dry delay:

[0144] ;

[0145] Among them, is the theoretical value of the tropospheric zenith dry delay, is the dry refractive index of each layer of the atmosphere;

[0146] The theoretical value of the tropospheric zenith wet delay is obtained by subtracting the theoretical value of the tropospheric zenith dry delay from the theoretical value of the tropospheric zenith total delay:

[0147] .

[0148] The method of forward modeling using the numerical integration method is as follows:

[0149] 1) Assume that the propagation path of the GNSS signal in the troposphere is a straight line, discretize the propagation path of the GNSS signal in the troposphere, divide the propagation path into several segments, and on the basis of the discretization of the propagation path, calculate the penetration points of each layer, that is, the points where the GNSS signal propagation path intersects perpendicularly with a specific height level. The position information of the penetration points is expressed as:

[0150] ;

[0151] Among them, is the position of the penetration point in the th layer, is the latitude coordinate of the penetration point, is the longitude coordinate of the penetration point, is the height of the penetration point;

[0152] 2) Calculate the atmospheric refractive index profile corresponding to each data point based on the constructed data points:

[0153] ;

[0154] Among them, represents the temperature of the th layer of the th data point, represents the water vapor pressure of the th layer of the th data point, represents the air pressure of the th layer of the th data point, is the atmospheric refractive index of the th layer of the th data point, then There are a total of layers, that is:

[0155] ;

[0156] 3) Accumulate the atmospheric refractive index along the signal propagation path through the path integral method to obtain the theoretical value of the tropospheric zenith total delay The calculation formula is:

[0157] ;

[0158] The theoretical value of the tropospheric zenith dry delay The calculation formula is:

[0159] ;

[0160] Among them, is the dry atmospheric refractive index, is the zenith dry delay above the upper limit height of the profile, Obtained from the Saastamoinen model:

[0161] ;

[0162] Among them, is the air pressure at the height of the layer penetration point, is the latitude at the height of the layer penetration point, is the height of the layer penetration point;

[0163] The theoretical value of the tropospheric zenith wet delay Is obtained by subtracting the theoretical value of the tropospheric zenith dry delay from the theoretical value of the tropospheric zenith total delay :

[0164] 。

[0165] In the above solution, the and calculation formulas are as follows:

[0166] ;

[0167] Among them, is the latitude of the GNSS receiver, is the longitude of the GNSS receiver, is the satellite azimuth corresponding to the th layer penetration point, The calculation formula of

[0168] ;

[0169] Among them, is the satellite elevation corresponding to the th layer penetration point, is the radius of the earth, is the height of the GNSS receiver, is the height of the th layer penetration point.

[0170] Step 4: Select one or more of the tropospheric zenith total delay, dry delay, and wet delay, use the optimization algorithm to invert the tropospheric atmospheric parameter profile, select the data point set with a smaller error between the theoretical value and the actual observation value for iteration, and finally obtain the optimal data points. After performing data quality control on the temperature, water vapor pressure, air pressure, and atmospheric refractive index corresponding to the optimal data points, output them as the final inversion result.

[0171] The optimization algorithm includes single-objective optimization or multi-objective optimization methods; among them, single-objective optimization uses genetic algorithms or particle swarm optimization algorithms, and the optimization objective is to minimize the root mean square error between the theoretical value and the actual observation value of one of the tropospheric zenith total delay, tropospheric zenith dry delay, and tropospheric zenith wet delay;

[0172] Multi-objective optimization uses non-dominated sorting genetic algorithms or multi-objective particle swarm algorithms, and its optimization objective is to minimize the root mean square error value 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 is expressed by the following formula:

[0173] ;

[0174] Among them, is the multi-objective function vector, is the root mean square error value between the actual observation value and the theoretical value of the zenith total delay, is the root mean square error value between the actual observation value and the theoretical value of the zenith dry delay, is the root mean square error value between the actual observation value and the calculated value of the zenith wet delay, represents the actual observation value of the total delay, represents the theoretical value of the total delay, represents the actual observation value of the dry delay, represents the theoretical value of the dry delay, represents the actual observation value of the wet delay; represents the theoretical value of the wet delay.

[0175] For a station with the co-location observation conditions of a ceilometer and a microwave radiometer, or the data reception conditions of a numerical weather prediction model output, the ceilometer, microwave radiometer or numerical weather prediction model output data are used as the constraint conditions for 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 (decreasing with height) above the cloud layer and a positive gradient (increasing with height) below the cloud layer; the temperature, pressure and water vapor pressure profiles output by the microwave radiometer or numerical weather prediction model are used as the initial values of the data constructed in the subsequent inversion algorithm;

[0176] If the station does not have the co-location observation conditions of a ceilometer and a microwave radiometer, and the data reception conditions of a numerical weather prediction model output, no constraint conditions are set.

[0177] When constraint conditions are set, the method for inverting the tropospheric atmospheric parameter profile using an optimization algorithm is as follows:

[0178] Use the constraint set to represent the external constraint condition set:

[0179] Use the cloud height data measured by the ceilometer to set the constraint conditions: the profile data are divided into a total of layers, and the height of each layer is , define the cloud height measured by the ceilometer corresponding to the profile layer number to obtain the layer number closest to the cloud height; according to the existing temperature and pressure to obtain the saturated water vapor pressure , then the constraint condition

[0180] is:

[0181] Among them, represents the th data point layer water vapor pressure, represents the th data point The saturated water vapor pressure of the layer is the layer corresponding to the cloud height indicating the number of layers corresponding to the cloud height of the water vapor pressure The number of layers corresponding to the cloud height of the saturated water vapor pressure;

[0182] Using the atmospheric parameter profiles output by the microwave radiometer or the numerical prediction model, select the external output profile that can be obtained in the actual situation from the two to set the constraint conditions; for the constructed data point profile and 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 not greater than a certain small value as the constraint condition, then the constraint condition is:

[0183] ;

[0184] Among them, represents the temperature of the constructed data point at the layer, represents the temperature of the external output profile at the layer, represents the water vapor pressure of the constructed data point at the layer, represents the water vapor pressure of the external output profile at the layer, represents the air pressure of the constructed data point at the layer, represents the air pressure of the external output profile at the layer, , and are the error thresholds of the allowable temperature, air pressure and water vapor pressure at the layer;

[0185] Then the final objective function is expressed as:

[0186] ;

[0187] Select the point set closest to the actual observation value from the data point set obtained by inversion, and through multiple rounds of iteration and gradual optimization, obtain the optimal data points. After performing data quality control on the temperature, water vapor pressure, air pressure and atmospheric refractive index corresponding to the optimal data points, output them as the inversion result.

[0188] The data quality control method is as follows:

[0189] For the case of multiple solutions, first eliminate the complex solutions; secondly, perform checks on climate threshold values, internal consistency of data, and spatial consistency. Based on the physical laws and statistical characteristics of the vertical distribution in the troposphere, eliminate the ambiguous solutions that do not conform to the objective physical environment, such as relative humidity exceeding 100%. Based on the logical relationships between meteorological elements, eliminate the solutions with contradictory relationships between meteorological elements. Based on the data of adjacent meteorological stations, eliminate the solutions with errors exceeding the reasonable range; obtain the optimal solution of the tropospheric atmospheric parameter profile for the final output.

[0190] The above description of the disclosed embodiments enables those skilled in the art to implement or use the present invention. Various modifications to these embodiments will be obvious to those skilled in the art, and the general principles defined herein can 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 these embodiments shown herein, but rather to the widest scope consistent with the principles and novel features disclosed herein.

Claims

1. A method for retrieving the tropospheric atmospheric parameter profile at sea based on GNSS, characterized in that, It includes the following steps: Step 1: Estimate the total zenith delay of the troposphere at the station location using GNSS real-time single-point positioning. Measure the meteorological elements at the station location using a meteorological element sensor. Calculate the dry zenith delay of the troposphere based on the meteorological elements. Separate the dry zenith delay of the troposphere from the total zenith delay of the troposphere to obtain the wet zenith delay of the troposphere, thereby obtaining the actual observed values of the total zenith delay, dry delay, and wet delay of the troposphere; Step 2: Perform cluster analysis using the historical sounding data of a sounding station near the station or the reanalysis data of the station location to obtain the statistical characteristics of the regional tropospheric atmospheric parameter profiles; Step 3: Construct data points regarding temperature, water vapor pressure, and air pressure using the statistical characteristics of the tropospheric atmospheric parameter profiles. Use the ray tracing method or numerical integration method combined with the atmospheric refractive index to perform forward calculations of the theoretical values of the total zenith 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 an optimization algorithm to invert the tropospheric atmospheric parameter profiles. Select a data point set with a smaller error between the theoretical value and the actual observed value for iteration. Finally, obtain the optimal data points. After performing data quality control on the temperature, water vapor pressure, air pressure, and atmospheric refractive index corresponding to the optimal data points, output them as the final inversion results; In Step 3, the method of performing forward calculation using the ray tracing method is as follows: First, taking the position of the GNSS antenna at the survey station as the starting position and the position vector , ray tracing is carried out in the zenith direction; Interpolate the vertical gradient of atmospheric refractive index according to the ray tracing calculation step size to obtain the atmospheric refractive index corresponding to each calculation step size ; Iteratively calculate from the th layer to the th layer until the cumulative height of the calculation step size reaches the upper limit of the inversion profile; Secondly, based on the layered atmosphere model, the delay of signal propagation is calculated. For each segment of the ray path, the propagation delay of the signal in this segment is calculated according to the refractive index and path length of this segment, and then the delays of each segment are accumulated to obtain the theoretical value of the total zenith delay at the upper limit of the inversion profile. Above the upper limit height of the profile, the water vapor is zero, and the total zenith delay is represented by the zenith dry delay through the Saastamoinen model. Then the theoretical value formula for the forward calculation of the total tropospheric propagation delay is as follows: where the water vapor is zero above the upper limit height of the profile, and the total zenith delay is represented by the zenith dry delay through the Saastamoinen model. Then the theoretical value formula for the forward calculation of the total tropospheric propagation delay is as follows: ； Among them, is the theoretical value of the tropospheric zenith total delay, is the atmospheric refractive index of each layer, is the ray at the layer path length, is the zenith dry delay above the upper limit height of the profile; ； Among them, is the air pressure at the upper limit height of the section, is the latitude at the upper limit height of the section, is the upper limit height of the section; Similarly, calculate the theoretical value of the dry zenith delay of the troposphere: ； Among them, is the theoretical value of the tropospheric zenith dry delay, is the dry refractivity of each layer of the atmosphere; The theoretical value of tropospheric zenith wet delay Obtained by subtracting the theoretical value of tropospheric zenith total delay from the theoretical value of tropospheric zenith dry delay: ； In Step 3, the method of performing forward calculation using the numerical integration method is as follows: 1) Assume that the propagation path of the GNSS signal in the troposphere is a straight line. Discretize the propagation path of the GNSS signal in the troposphere, divide the propagation path into several segments. Based on the discretization of the propagation path, calculate the penetration points of each layer, that is, the points where the GNSS signal propagation path intersects perpendicularly with a specific altitude layer. The position information of the penetration points is expressed as: ； Among them, is the position of the layer-piercing point of the layer, is the latitude coordinate of the layer-piercing point, is the longitude coordinate of the layer-piercing point, is the height of the layer-piercing point; 2) According to the constructed data points, calculate the atmospheric refractive index profile corresponding to each data point: ； Among them, represents the temperature of the th data point at the th layer, represents the vapor pressure of the th data point at the th layer, represents the air pressure of the th data point at the th layer, is the atmospheric refractive index of the th data point at the th layer, then There are layers, that is: ； 3) By using the path integral method to accumulate the atmospheric refractive index along the signal propagation path, the theoretical value of the total tropospheric zenith delay is obtained. The calculation formula is as follows: ； The theoretical value of tropospheric zenith dry delay The calculation formula is as follows: ； wherein, is the atmospheric dry refractivity, is the zenith dry delay above the upper height of the profile, obtained from the Saastamoinen model: ； Among them, is the air pressure at the height of the layer-piercing point , is the latitude at the height of the layer-piercing point , is the height of the layer-piercing point; Theoretical value of tropospheric zenith wet delay From the theoretical value of the total tropospheric zenith delay And the theoretical value of the tropospheric zenith dry delay By taking the difference: 。 2. The method for inverting the tropospheric atmospheric parameter profile based on GNSS according to claim 1, wherein For a station with the co-location observation conditions of a ceilometer and a microwave radiometer, or the data reception conditions of a numerical weather prediction model output, use the ceilometer, microwave radiometer, or numerical weather prediction model output data as the constraint conditions for 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 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 weather prediction model are used as the initial values of the data constructed in the subsequent inversion algorithm; If the station does not have the co-location observation conditions of a ceilometer and a microwave radiometer, as well as the data reception conditions of a numerical weather prediction model output, then no constraint conditions are set.

3. The method for inverting the tropospheric atmospheric parameter profile based on GNSS according to claim 1, characterized in that, In Step 1, the method of obtaining the actual observed values of the total zenith delay, dry delay, and wet delay of the troposphere is as follows: (1) The dual-frequency raw observables are obtained by using a GNSS receiver. Combining with the Beidou satellite-based broadcast PPP-B2b corrections, the ionospheric delay error is eliminated through the ionosphere-free combination, and the actual observation value of the tropospheric zenith total delay is extracted. ; (2) Calculate the dry zenith delay of the troposphere using meteorological elements: ； Among them, is the actual observed value of the tropospheric zenith dry delay, is the geographical latitude of the station, is the air pressure at the station location, is the station elevation; (3) Calculate the wet zenith delay of the troposphere using the total zenith delay of the troposphere and the dry zenith delay of the troposphere: ； Among them, is the actual observed value of the tropospheric zenith wet delay.

4. A method for inverting the tropospheric atmospheric parameter profile based on GNSS according to claim 1, characterized in that The method adopted for the cluster analysis is the 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 profiles is as follows: Select the historical radiosonde data of the radiosonde station near the measurement station or the reanalysis data of the measurement station location to construct a historical profile dataset; preprocess the temperature, relative humidity, and refractive index profiles in the dataset: first, detect and remove outliers from the original profile data, and then calculate the value range of meteorological parameter data, the first-order change rate and the second-order change rate of the parameter with height for each vertical layer of each profile. Normalize the data to unify the parameters and their change rates in different height layers. Group the normalized profile dataset by month, perform nonlinear dimensionality reduction and topological mapping on the dataset, train the SOM self-organizing mapping network, and iteratively adjust through competitive learning and neighborhood functions. For each input data, calculate its Euclidean distance from each neuron node, cluster the data with similar characteristics into the same neuron node. After meeting the iteration termination condition, the training is completed, and each data is assigned to the cluster represented by the neuron most similar to it to obtain the final clustering result; calculate the statistical characteristics for each cluster. First, extract the cluster center, and then statistically calculate the value range of meteorological parameters layer by layer along the height dimension. Further extract the first-order change rate and the second-order change rate of each meteorological element to obtain the in-class statistical distribution characteristics, and finally obtain the statistical characteristics of the tropospheric atmospheric parameter profile.

5. A method for inverting the tropospheric atmospheric parameter profile based on GNSS according to claim 1, characterized in that In step 3, the specific method for constructing data points regarding temperature, water vapor pressure, and air pressure using the statistical characteristics of the tropospheric atmospheric parameter profile is as follows: The assumed constructed data points are , where represents temperature,[[]] represents water vapor pressure,[[]] represents air pressure,[[]] represents the th data point,[[]] represents the temperature profile,[[]] represents the water vapor pressure profile,[[]] represents the air pressure profile; For the temperature profile , is the total number of layers in the profile, denotes the -th layer in the layer data; assuming that the minimum value of a single layer in the profile is , the maximum value is , the first-order change rate of the profile limited between the minimum and maximum values due to the profile characteristics is , and its value range is defined between the minimum value and the maximum value ; the second-order change rate is , and its value range is defined between the minimum value and the maximum value ; For the first-layer random number , its generation formula is: ； Here, represents a uniformly distributed random variable on the interval to ensure that the value is within the interval ; Subsequent random number , for the first-order change rate limit, by defining the first-order change rate of each layer , the following is obtained by transformation: ； Among them, is the height difference between the upper and lower layers of data; For the second-order rate-of-change limit without considering the height variable, it is defined as , and by transformation, we get: ； Take the intersection of the first-order and second-order change rate limit conditions for the random numbers in the above two formulas, and finally obtain the value range as follows: ； If this intersection is empty, then take the random number generated for the first-order rate-of-change limit, i.e., , denotes a uniformly distributed random variable on the interval ; According to the above process, based on the statistical characteristics of the single-layer maximum and minimum values of the given profile and the random change rate, a temperature profile that meets the requirements is generated. A random number sequence is generated, and similarly, a water vapor pressure profile sequence and an air pressure profile sequence are generated to form which are used as data points.

6. The method for inverting the tropospheric atmospheric parameter profile based on GNSS according to claim 1, characterized in that For each penetration point and The calculation formula is as follows: ； Wherein, is the latitude of the GNSS receiver, is the longitude of the GNSS receiver, is the satellite azimuth corresponding to the layer-piercing point of the layer, and its calculation formula is as follows: ； Among them, is the satellite elevation angle corresponding to the layer penetration point, is the radius of the earth, is the height of the GNSS receiver, is the height of the layer penetration point.

7. A method for retrieving the tropospheric atmospheric parameter profile based on GNSS according to claim 1, characterized in that, In step 4, the optimization algorithm includes single-objective optimization or multi-objective optimization methods; among them, single-objective optimization uses genetic algorithms or particle swarm optimization algorithms, 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, tropospheric zenith dry delay, and tropospheric zenith wet delay. Multi-objective optimization uses non-dominated sorting genetic algorithms or multi-objective particle swarm algorithms, and its optimization goal is to minimize the root mean square error values between the theoretical values and the actual observed values of at least two of the tropospheric zenith total delay, tropospheric zenith dry delay, and tropospheric zenith wet delay, which is expressed by the following formula: ； Among them, is the multi-objective function vector, is the root mean square error value between the actual observation value and the theoretical value of the total zenith delay, is the root mean square error value between the actual observation value and the theoretical value of the dry zenith delay, is the root mean square error value between the actual observation value and the calculated value of the wet zenith delay, represents the actual observation value of the total delay, represents the theoretical value of the total delay, represents the actual observation value of the dry delay, represents the theoretical value of the dry delay, represents the actual observation value of the wet delay; represents the theoretical value of the wet delay.

8. A method for retrieving the tropospheric atmospheric parameter profile based on GNSS according to claim 2, characterized in that, In step 4, when setting constraint conditions, the method for inversely calculating the tropospheric atmospheric parameter profile using the optimization algorithm is as follows: Use a set of constraints to represent an external set of constraint conditions: Set constraint conditions using cloud height data measured by a ceilometer: The profile data is divided into a total of layers, and the height of each layer is . Define the profile layer number corresponding to the cloud height measured by the ceilometer to obtain the layer number closest to the cloud height; Obtain the saturated water vapor pressure from the existing temperature and air pressure . Then the constraint condition is: ； Among them, represents the water vapor pressure of the th data point at the th layer, represents the saturated water vapor pressure of the th data point at the th layer, is the layer corresponding to the cloud height, represents the water vapor pressure of the layer number corresponding to the cloud height , represents the saturated water vapor pressure of the layer number corresponding to the cloud height ; Using the atmospheric parameter profiles output by a microwave radiometer or a numerical weather prediction model, select the external output profiles that can be obtained in the actual situation from the two to set the constraint conditions; for the constructed data point profiles and the output profiles , set the absolute value of the difference in temperature and water vapor pressure data for each layer of the two to be not greater than a certain small value as the constraint condition, then the constraint condition is: ； Among them, represents the temperature of the structural data point at the layer, represents the temperature of the external output profile at the layer, represents the water vapor pressure of the structural data point at the layer, represents the water vapor pressure of the external output profile at the layer, represents the air pressure of the structural data point at the layer, represents the air pressure of the external output profile at the layer, , and are the error thresholds of the allowable temperature, air pressure and water vapor pressure at the layer; Then the final objective function is expressed as: ； Select the point set closest to the actual observed value from the inversely calculated data point set, and through multiple rounds of iteration and gradual optimization, obtain the optimal data point. After performing data quality control on the temperature, water vapor pressure, air pressure, and atmospheric refractive index corresponding to the optimal data point, output them as the inversion result.

Images