Beidou tomography of tropospheric water vapor based on lstm and ray ordered subset constraint strategy

Abstract

The application discloses a Beidou troposphere water vapor tomography method based on an LSTM and a ray ordered subset constraint strategy, first constructs a water vapor density prediction model, and uses the water vapor density prediction model to predict the water vapor density distribution of a tomography target time as an iteration initial field; then, based on Beidou observation data, an effective ray set and corresponding slant path water vapor values of the tomography target time are obtained, the intercepts of each ray in a preset three-dimensional voxel grid are calculated to construct an intercept matrix, and all effective rays are divided into multiple ordered balanced subsets, so as to construct an effective ray library; finally, a fusion algebra reconstruction algorithm is adopted, and the effective ray library is used for iteration update, a constructed Beidou troposphere tomography matrix equation is solved, and a three-dimensional water vapor density inversion result is output. The method can effectively overcome problems such as insufficient tomography inversion accuracy, low inversion calculation efficiency and weak robustness.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of navigation satellite remote sensing inversion, and particularly relates to a Beidou troposphere water vapor tomography method based on LSTM and ray ordered subset constraint strategy. BACKGROUND

[0002] The troposphere is the low layer of atmosphere from the ground to about 8 to 18 km, is the most concentrated area of atmospheric quality and water vapor, and is also the occurrence layer of main weather systems. The water vapor and temperature thereof are controlled by the ground heating and circulation, and present spatiotemporal non-uniformity. The GNSS (Global Navigation Satellite System) signal propagates in the troposphere to produce refraction and delay, especially the wet delay which is difficult to model and changes dramatically, and is an important error source of high-precision positioning. The troposphere water vapor is the main energy and indicator factor of weather evolution, especially strong convective weather (heavy rain, strong convective storm, etc.), has strong spatiotemporal non-uniformity, and is crucial to weather forecast, climate research, and disaster prevention and reduction. Therefore, the construction of troposphere three-dimensional (even four-dimensional) high spatiotemporal resolution and high-precision water vapor structure by using GNSS and other multi-source observations not only helps to reduce satellite navigation positioning errors and improve precise measurement and timing accuracy, but also provides fine water vapor information support for the monitoring and early warning of small-scale disastrous weather, and is an important problem to be solved in current troposphere water vapor tomography.

[0003] The existing troposphere tomography methods can be roughly divided into four categories:

[0004] The first category is to use radio sounding and ground-based profiler instruments. The sounding balloon provides temperature, humidity, and pressure vertical profile, and the ground-based microwave radiometer and laser radar continuously inverts water vapor or aerosol profile at fixed sites, which is equivalent to providing one-dimensional vertical information at limited sites and reconstructing three-dimensional structure by interpolation, and is the main information source of early three-dimensional water vapor analysis and tomography. This kind of observation is direct and high in precision, but the cost is high, and it is difficult to meet the demand of local high spatiotemporal resolution three-dimensional tomography.

[0005] The second category of methods uses satellite-borne infrared, microwave radiometer and airborne passive remote sensing to obtain temperature, humidity and cloud water content, and combines radiation transfer to invert the three-dimensional field of the troposphere; GNSS occultation and atmospheric refraction tomography also belong to this category. Multi-orbit, multi-angle observation has tomography geometry, and can construct three-dimensional refractive index / water vapor field on a global or regional scale, has wide coverage and good representation of ocean structure, but is limited by orbit and observation geometry and inversion model, and has limited representation of local structure and low spatial resolution.

[0006] The third method is to obtain continuous four-dimensional atmospheric analysis field by assimilating sounding, ground-based remote sensing, satellite, radar and other multi-source observations into the model through variational assimilation or ensemble Kalman filter. This method has good spatio-temporal continuity and strong physical consistency, but it is highly dependent on the model and assimilation scheme, and the local small-scale water vapor structure is still prone to smoothing.

[0007] The fourth method is to use ground-based GNSS signal observation data to develop tropospheric tomography. The troposphere is discretized into voxel or node grid, and the tomography equation is established based on GNSS slant path wet delay or slant path water vapor value. The result of slant path water vapor value inversion is water vapor density, and the result of slant path wet delay inversion is wet refractive index. Combined with the iterative initial value, the equation is solved by algebraic reconstruction, Kalman filter and other methods to obtain real-time three-dimensional water vapor structure of the troposphere.

[0008] Compared with other methods, ground-based GNSS tomography has low cost, all-weather and all-day, high time resolution, and can make up for the shortcomings of sounding and satellite in near-surface layer and local monitoring. However, the existing methods for obtaining the initial value of iterative inversion of water vapor tomography are all based on experience, fitting water vapor density model, radiosonde data or general water vapor density prediction model, which has the problems of poor estimation of water vapor distribution in severe weather, low time resolution, poor real-time performance, delay in sensing water vapor changes in extreme weather, etc., which cannot obtain accurate water vapor density as the initial value of tomographic iteration, and further affect the accuracy of Beidou tropospheric tomography results. In addition, the commonly used iterative method for solving tomographic equation needs to process all effective rays in each iteration, which has the problems of low calculation efficiency, slow convergence of iteration, large amount of calculation, etc. The commonly used iterative algorithms such as Multiplicative Algebraic Reconstruction Technique (MART) and Simultaneous Iterative Reconstruction Technique (SIRT) are prone to overfitting of observation rays, resulting in excessive updating, noise amplification and iterative shock. MART has too many iterations, which can amplify noise, while SIRT has too strong overall smoothness, resulting in loss of details. At the same time, the convergence efficiency and result stability of this kind of method are limited, the same and fixed relaxation factor is used for all voxel layers, which cannot adapt to the different degrees of ill-conditioning of high and low voxel layers, and iterative shock occurs, and even when the top layer is constrained, the error will be propagated to the middle and lower layers of water vapor.

[0009] Therefore, researching a tomographic water vapor initial value model with good real-time performance, high spatiotemporal resolution, and strong water vapor sensing capability, along with a suitable ray processing strategy and an algebraic reconstruction method that features fast iterative convergence, high computational efficiency, adaptability to different ill-conditioning levels in different voxel layers, and the advantages of both MART and SIRT, can effectively alleviate the problems of lack of accurate and reliable initial value models, excessive inversion iterations, and insufficient inversion accuracy in BeiDou tropospheric tomography. This will significantly improve inversion accuracy and reliability, while also exhibiting good real-time performance and robustness. Consequently, it will provide more robust and efficient technical support for tropospheric water vapor tomography. Summary of the Invention

[0010] The purpose of this invention is to overcome the shortcomings of the prior art and propose a BeiDou tropospheric water vapor tomography method based on LSTM (Long Short Term Memory) and ray ordered subset constraint strategy. This method can effectively overcome the problems of insufficient tomographic inversion accuracy, low inversion calculation efficiency, and weak robustness caused by the lack of accurate and reliable initial value models, slow iteration convergence, large amount of computation, iteration oscillation, and excessive reverse update.

[0011] To achieve the above objectives, the technical solution specifically adopted by the present invention is as follows:

[0012] A BeiDou-based tropospheric water vapor tomography method based on LSTM and ray ordered subset constraint strategy includes the following steps:

[0013] Step 1: Construct a water vapor density prediction model, and use the water vapor density prediction model to predict the water vapor density distribution at the target time of the chromatography as the initial field for iteration;

[0014] Step 2: Based on BeiDou observation data, obtain the effective ray set and its corresponding oblique path water vapor value at the time of the tomographic target. Calculate the intercept of each ray in the preset three-dimensional voxel grid to construct the intercept matrix. Divide all effective rays into multiple ordered balanced subsets to construct an effective ray library containing the oblique path water vapor value, intercept matrix and ordered balanced subset of effective rays, and save it in text form.

[0015] Step 3: Using a fusion algebraic reconstruction algorithm, combined with dual-source relaxation factor constraints based on height prior information and iterative residual information, and using the effective ray library for iterative updates, the BeiDou tropospheric tomography matrix equation constructed with the iterative initial field, intercept matrix and oblique path water vapor value is solved, and the three-dimensional water vapor density inversion result is output.

[0016] Preferably, the steps for constructing the water vapor density prediction model include:

[0017] Obtain radio sounding water vapor data of historical period, and calculate water vapor density and atmospheric precipitable water vapor (PWV) of multiple height layers at multiple historical time points based on the radio sounding water vapor data;

[0018] Construct a neural network model comprising an input feature layer, an LSTM sequence modeling layer, a ReLU activation function layer, a full connection layer, and a loss regression layer;

[0019] Train the neural network model by taking a combination of water vapor density sequence, atmospheric precipitable water vapor (PWV), and time encoding features of historical time points as input and taking water vapor density of a target time point as output, to obtain the water vapor density prediction model.

[0020] Preferably, the time encoding feature is a feature value obtained by sine function and cosine function transformation on the accumulated day of the target time point.

[0021] The step of predicting the water vapor density distribution of the target time point as the initial field of iteration comprises: inputting the water vapor density sequence of N W time steps before the tomography time point, the atmospheric precipitable water vapor (PWV) at the tomography time point (obtained in step 2), and the tomography time point encoding feature into the water vapor density prediction model for prediction, and taking the predicted water vapor density sequence as the initial field of iteration.

[0022] Preferably, in step 2, the step of obtaining the effective ray set and the slant path water vapor value comprises:

[0023] Calculate zenith total delay (ZTD) of each receiver site based on Beidou observation data, and calculate zenith hydrostatic delay (ZHD) and zenith wet delay (ZWD) in combination with site meteorological parameters;

[0024] Convert the zenith wet delay (ZWD) into atmospheric precipitable water vapor (PWV) by using a conversion factor, and project the zenith wet delay (ZWD) to each slant path by using a mapping function and a horizontal gradient model in combination with the elevation and azimuth of each ray, to obtain the slant path water vapor value (SWV) of each ray, and further construct an SWV matrix.

[0025] Preferably, the construction method of the intercept matrix is:

[0026] First, convert the geometric relationship between the receiver and the satellite into a three-dimensional voxel grid in the local rectangular coordinate system, express each signal path by using a ray parameter equation, and calculate the intersection points of each voxel boundary, calculate the straight line segment length of the ray in the voxel according to the start point and end point parameters of the ray in the single voxel, and the length is the intercept corresponding to the voxel; and further construct the intercept matrix.

[0027] Preferably, the step of dividing the effective rays into multiple ordered balanced subsets comprises:

[0028] The effective rays are evenly distributed by mixing according to at least one of the elevation angle, the azimuth angle, the receiver station to which the ray belongs, and the path length or the number of voxels crossed by the ray, so that each subset contains rays with different attribute ranges and each subset is representative of the statistical characteristics of the entire ray set.

[0029] Preferably, the fusion algebraic reconstruction algorithm comprises a periodic alternating iterative algorithm of the multiplicative algebraic reconstruction technique (MART) and the simultaneous iterative reconstruction technique (SIRT).

[0030] Preferably, the Beidou tropospheric water vapor tomography method based on the LSTM and the ordered subset constraint strategy of the ray comprises a periodic alternating iterative manner, that is, after completing a preset number of MART iterations on all ordered balanced subsets, performing one SIRT iteration based on all rays.

[0031] Preferably, the relaxation factor constraint based on the high-priority information is:

[0032] Different relaxation factor weights are given to voxels at different altitudes, wherein for voxels with an altitude higher than a preset threshold, the weight of the voxel is attenuated in a negative exponential relationship according to the altitude difference exceeding the threshold; and for voxels with an altitude lower than or equal to the preset threshold, the weight of the voxel is a reference value.

[0033] Preferably, the relaxation factor constraint based on the iterative residual information is:

[0034] After each iteration update, the residual between the current inversion obtained oblique path water vapor value and the actual oblique path water vapor value is calculated, and the relaxation factor size of the next iteration is dynamically adjusted according to the statistical quantity of the residual.

[0035] Preferably, the way of dynamically adjusting the relaxation factor according to the statistical quantity of the residual is: inputting the average absolute residual value of the residual into a hyperbolic tangent function for mapping, and taking the mapping result as an adjustment coefficient to adjust the relaxation factor.

[0036] Preferably, after solving the tomography equation iteratively, an empirical constraint in the horizontal and vertical directions is introduced to constrain the water vapor density value in voxels not crossed by rays.

[0037] Preferably, the rays with an elevation angle greater than a set angle threshold and a projection height in the zenith direction greater than a set height threshold are selected as the effective rays.

[0038] The application mainly provides a Beidou tropospheric water vapor tomography method based on an LSTM fusion PWV water vapor density prediction model, a ray ordered subset constraint strategy and a double-source relaxation factor constrained fusion algebraic reconstruction method (DSCR-HART). The water vapor density prediction model based on the LSTM fusion PWV can provide accurate, reliable and high space-time resolution initial values for tomography for iterative inversion, has good sensing ability for water vapor changes in bad weather and can significantly improve the tomography inversion accuracy. The ray subset constraint strategy divides the satellite rays obtained in the sampling window into multiple subsets, each subset has overall representativeness for iteration, and the information of multiple ordered balanced subsets is combined to effectively reduce the calculation amount, improve the iteration convergence rate and effectively alleviate the problem of large calculation amount caused by traversing all rays in each iteration. The double-source constraint on the relaxation factor can effectively adapt to the ill-condition degree of different height layers of voxels, effectively improve the sensing ability of water vapor vertical motion distribution, ensure fast convergence speed, suppress noise, prevent reverse update and reduce inversion error. The MART and SIRT fusion method has the advantages of both methods, has stronger anti-noise ability, better robustness, appropriate smoothing performance and can timely correct the noise voxels that may be introduced by the smoothing MART, and has high iteration efficiency. This fusion iteration mode cooperates with the double-source constrained relaxation factor, the ordered balanced subset division constraint strategy and the excellent water vapor density prediction model to provide initial values for tomography iterative inversion, and can further improve the reliability and accuracy of the Beidou tropospheric tomography.

[0039] The application has the following characteristics and beneficial effects:

[0040] The method described in this invention can fully utilize BeiDou Navigation Satellite System data and radiosonde data to perform accurate water vapor tomography of the troposphere. Compared with existing GNSS tropospheric tomography methods based on empirical water vapor density models as iterative initial values, traditional tomography based on radiosonde data, and GNSS occultation methods, it solves the problems of inaccurate empirical models, incorrect estimation of vertical water vapor distribution leading to large tomography errors in severe weather, the lack of horizontal resolution and extremely low temporal resolution of radiosonde data, and the low resolution of GNSS occultation methods. Compared with GNSS tropospheric tomography methods based on traditional algebraic reconstruction methods, it solves the problems of low computational efficiency and slow iterative convergence of traditional methods, exhibiting excellent convergence rate, low computational load, and good stability and accuracy. The water vapor density prediction model fused with LSTM and PWV fully leverages the advantages of radiosonde data and BeiDou data, enhancing the ability to sense water vapor movement. This effectively addresses the problems of large errors in empirically fitted models, poor water vapor sensing capabilities in low-severity weather, and delayed sensing in general prediction models when encountering sudden changes. It obtains high-resolution and accurate initial values ​​for tomographic water vapor density iterations, significantly improving tomographic performance. By dividing rays into multiple balanced subsets sequentially, each subset can represent the whole and iterate sequentially, suppressing excessive back-updates, resulting in less computation and faster, more accurate convergence per generation. The dual-source constraint relaxation factor effectively adapts to different levels of ill-conditioning in voxel layers at high altitudes with low water vapor density, forming a progressive soft constraint in the vertical direction. This better suppresses instability caused by ill-conditioning in high-level data while preventing the propagation of top-level forced constraint errors to mid- and low-level water vapor-dense areas. Dynamically adjusting the iterative relaxation factor based on SWV residuals ensures rapid convergence while effectively suppressing noise, resulting in smaller water vapor density inversion errors and stronger sensing capabilities for vertical water vapor movement. This method addresses the challenges of adjusting fixed relaxation factors and adapting to varying degrees of ill-conditioning at different altitudes and voxels. By employing a fusion iterative approach combining MART and SIRT—using MART for fast approximation and SIRT for periodic stability—it leverages the strengths of both algorithms to balance local and global iterative adjustments. This reduces noisy voxels, avoids excessive smoothing, and effectively mitigates the rank deficiency problem during inversion, resulting in more stable BeiDou tropospheric tomography results. This further ensures the accuracy, real-time performance, robustness, and reliability of BeiDou tropospheric water vapor tomography. Attached Figure Description

[0041] Figure 1 This is a flowchart illustrating the process of acquiring water vapor data and constructing a water vapor density prediction model based on LSTM fusion PWV in an embodiment of the present invention.

[0042] Figure 2A flow chart for constructing an effective ray library and an intercept matrix based on Beidou observation data in the embodiment of the present application is provided.

[0043] Figure 3 A flow chart for predicting the water vapor density at the tomographic time point as an initial value of iteration by a prediction model in the embodiment of the present application is provided.

[0044] Figure 4 A flow chart for solving the tomographic equation to output the water vapor density based on the ray ordered subset constraint strategy and the DSCR-HART in the embodiment of the present application is provided. DETAILED DESCRIPTION

[0045] The present application will be described in detail below with specific embodiments. The following embodiments will help those skilled in the art to further understand the present application, but do not limit the present application in any form. It should be noted that the embodiments in the present application and the features in the embodiments can be combined with each other without conflict.

[0046] The application aims to provide a Beidou tropospheric water vapor tomography method based on an LSTM fusion PWV water vapor density prediction model, a ray ordered subset constraint strategy and a double-source relaxation factor constrained fusion algebraic reconstruction method, thereby effectively improving the accuracy, real-time performance, stability and reliability of Beidou tropospheric water vapor tomography. The LSTM tomographic water vapor density prediction model fused with real-time Beidou PWV has strong ability to perceive water vapor content and vertical water vapor distribution in severe weather, so that the iteration initial value of the tomographic equation is accurate and reliable, and the tropospheric water vapor structure can be correctly inverted. The strategy of dividing the ray ordered balanced subset constraint further optimizes the calculation process, reduces the calculation complexity by merging the information of multiple balanced subsets, and improves the iteration convergence rate. The double-source constraint iteration relaxation factor can cope with the degree of ill-conditioning of different voxel layers. For the voxel layer with high water vapor content in the middle and low layers, no restriction is made, while for the voxel layer with generally low water vapor content in the high layer, the iteration change step is reduced to form a kind of soft constraint, which can avoid the ill-conditioning of the high layer and the downward propagation of constraint errors. At the same time, the iteration rate can be adjusted adaptively according to the residual matrix, thereby improving the convergence speed and accuracy of inversion and avoiding over-inversion, and the accuracy of reconstructed tropospheric water vapor structure is improved. HART can effectively suppress noise and have moderate smoothness by introducing periodic SIRT based on MART, and the iteration is fast and stable, which can effectively overcome the problems of excessive fitting of observation rays, excessive update, noise amplification and iteration shock. The LSTM fusion PWV water vapor prediction model can provide the tomographic equation with real-time, high spatial and temporal resolution, and stable and reliable water vapor density iteration initial value. The ray ordered subset constraint strategy can effectively reduce the calculation amount and reduce the calculation complexity. DSCR-HART has significant advantages in improving inversion accuracy, accelerating convergence and improving calculation efficiency. The fusion of the three can better cope with the challenges in tropospheric tomography. Therefore, the application solves the problems of insufficient accuracy of the existing GNSS tropospheric tomographic initial value model, poor real-time performance, delayed perception, weak vertical water vapor motion perception ability and the like, and effectively alleviates the problems of large overall calculation amount and large number of iterations. At the same time, it effectively overcomes the problems of low tomographic inversion result accuracy, poor robustness and reliability caused by the fixed relaxation factor, overall ray update, single iteration update mode and excessive update of the existing traditional method. The application can more accurately describe the complex atmospheric water vapor distribution by optimizing the tomographic water vapor density iteration initial value model, the ray subset constraint strategy and the optimized iteration algorithm, thereby improving the application performance of Beidou water vapor tomography in high-precision positioning, weather forecasting and disaster monitoring. It can effectively reduce the water vapor delay error and improve the GNSS positioning accuracy. In addition, the enhanced vertical water vapor motion perception ability improves the accuracy of weather forecasting and provides more reliable decision support for the early warning and emergency response of natural disasters such as rainstorms and floods, comprising the following steps:

[0047] A, a water vapor density prediction model is constructed, and the water vapor density distribution of the tomographic target time is predicted as an initial field for iteration using the water vapor density prediction model. As shown in Figure 1 、 Figure 3 , the steps include step 1 and step 2:

[0048] Step 1: Obtain historical radiosonde water vapor data and divide the voxels of the tomographic region.

[0049] First, the historical sounding data provided by the radiosonde data website is downloaded, and the water vapor density at different altitudes in accordance with the voxel division strategy is calculated and processed according to the sounding data combined with the corresponding formula. At the same time, the PWV is calculated according to the formula, and the specific formula is as follows:

[0050]

[0051]

[0052] In the formula, is the water vapor density , is the water vapor mixing ratio , is the atmospheric pressure , is the specific gas constant of water vapor, and is taken as , is the temperature . is the atmospheric precipitable water , is the local gravitational acceleration , is the total number of atmospheric layers in the sounding data, represents the th atmospheric layer in the sounding data and are the pressures of the th and the th layers, respectively and and are the specific humidities of the th and the th layers, respectively.

[0053] In this embodiment, a certain region is selected as the chromatography experiment area, and the voxel division strategy is: east longitude 113.88°-114.36°, interval 0.06°; north latitude 22.18°-22.6°, interval 0.06°; the altitude is unevenly divided, 0 km-0.5 km-1 km-1.7 km-2.4 km-3.1 km-4 km-4.9 km-5.8 km-6.9 km-8 km-9.1 km-10.2 km. The certain region is divided into an 8x7x12 three-dimensional grid, and the longitude and latitude are mapped to the distance through the local plane approximation model of the WGS-84 ellipsoid and the radius of the earth to form a local rectangular coordinate system. The historical data of the 45004 sounding station from 2018 to 2024 are selected for processing, and the data are standardized by Z-Score standardization, and then input into the LSTM network algorithm constructed in the next step for prediction training.

[0054] Step 2: Constructing the water vapor density prediction model of LSTM fusion PWV and other parameters and completing the training of the model.

[0055] The model in this embodiment is mainly composed of five parts in order: input feature layer, LSTM sequence modeling layer (hidden unit number is 15), ReLU activation layer, full connection layer (12-dimensional output) and loss regression layer. The input feature layer includes time coding feature value , , PWV and 12 height point (obtained in step 1) water vapor density sequence, a total of 15-dimensional input features, is the annual accumulated day of the target prediction time point, is the total number of days in a year, and the time sin-cos mapping is used to eliminate the discontinuity of the annual accumulated day value and explicitly encode the annual cycle phase information, which is conducive to the network to stably and fully learn the seasonal variation characteristics, and PWV is the PWV of the target prediction time point. The LSTM layer expands the input historical feature sequence in time and models the gate memory, extracts key time sequence information, and realizes the selection and update of information through input gate, forget gate and output gate, etc. The activation function of each gate is , and the activation function of the candidate cell state is , which all belong to the common sense in the field of neural networks and will not be described in detail. The ReLU activation layer performs element-wise nonlinear transformation on the output of the LSTM, thereby enhancing the network's ability to fit complex nonlinear relationships. The full connection layer maps the high-dimensional features obtained by the ReLU layer to a 12-dimensional output space through linear transformation, realizing the weighted mapping from time sequence features to each output variable. The loss regression layer constructs a loss function according to the error between the network prediction value and the true value, which is used to guide the back propagation and update of the parameters, thereby completing the regression learning of multiple output continuous variables.

[0056] The model adopts LSTM to capture the input sequence dependency, ReLU to introduce nonlinearity, enhance the network's fitting ability for complex features, and define a fully connected layer to map the input 15-dimensional feature data output to the output target 12-dimensional. According to the size of the data, select the training principle of small batch and multiple times, in this case, the number of samples used for each training iteration is 50, and the network is trained. At the same time, select Adam (Adaptive Moment Estimation, adaptive moment estimation) as the solution optimizer, which combines the advantages of Momentum (momentum) and RMSProp (root mean square propagation), has fast convergence speed and is not sensitive to learning rate, and is the mainstream optimizer of deep learning at present. Gradient clipping is selected to limit the maximum value of the gradient to prevent gradient explosion.

[0057] The data obtained in step 1 is input into the network for training, and finally the water vapor density prediction model based on LSTM fusion PWV is built. For the accuracy of the model, in this embodiment, the daily average RMSE (Root Mean Square Error, root mean square error) and MAE (Mean Absolute Error, mean absolute error) of the water vapor density from January to December 2025 obtained by different models are used to represent.

[0058]

[0059] Table 1

[0060] In Table 1, BP represents fitting historical water vapor density using BPNN (Back Propagation Neural Network, error back propagation neural network) to obtain water vapor density data from January to December 2025; BP-PWV represents fitting historical water vapor density using BPNN fusion historical PWV to obtain water vapor density data from January to December 2025; LSTM represents predicting water vapor density using LSTM combined with historical water vapor density data; and LSTM-PWV represents predicting water vapor density using LSTM fusion real-time PWV combined with historical water vapor density data.

[0061] As can be seen from Table 1, the water vapor density prediction model based on LSTM fusion PWV has the highest accuracy and the best effect, and the average RMSE of the other three fitting and general prediction models is improved by 43.71%, 21.01%, and 12.15%, respectively, and the MAE is improved by 43.80%, 26.47%, and 12.82%, respectively.

[0062] B. Based on the Beidou observation data, the effective ray set at the tomographic target time and the corresponding slant path water vapor value are obtained, the intercept of each ray in the preset three-dimensional voxel grid is calculated to construct an intercept matrix, and all effective rays are divided into multiple ordered balanced subsets, thereby constructing an effective ray library containing slant path water vapor values of effective rays, intercept matrices and ordered balanced subsets, and saving in text form. As shown in Figure 2 , specifically comprising steps 3 and 4:

[0063] Step 3: Calculate the SWV of the ray according to the Beidou observation data at the tomographic time point. Collect the Beidou original observation data at the tomographic time point, extract the ZTD of all receiving stations in the tomographic area through the Beidou data, there are 18 receiving stations in the tomographic area, combine the geographical positions of different receiving stations to calculate the ZHD according to the Sa model, the Sa model is an empirical model algorithm proposed by Saastamoinen to calculate ZHD. Subtract ZHD from ZTD to finally solve ZWD, the calculation formula is as follows:

[0064]

[0065]

[0066] In the formula, is the zenith hydrostatic delay , is the zenith hydrodynamic delay , is the zenith total delay , is the atmospheric pressure at the Beidou receiver , is the altitude of the Beidou receiver station . is the latitude of the Beidou receiver station.

[0067] ZWD and the conversion factor can be used to calculate SWV, and the conversion factor calculation formula is as follows:

[0068]

[0069] In the formula, is the conversion factor, is the liquid water density constant, taking , is the water vapor specific gas constant, taking , and are atmospheric refraction constants, taking and , is the atmospheric weighted mean temperature The empirical algorithm can be used to calculate, which is common knowledge in the field, The meanings of the formulas are the same.

[0070] To calculate SWV, different elevation and azimuth angles of the ray are combined, the mapping function is combined to compensate for the difference in the azimuth of the atmosphere, and the calculation formula is as follows:

[0071]

[0072]

[0073]

[0074] In the formula, is the slant path wet delay , is the slant path water vapor value , is the zenith wet delay, is the conversion factor, is the mapping function, and the mapping function is selected is the north-south horizontal gradient, is the east-west horizontal gradient, The mapping function and the gradient both belong to common knowledge in the field, and will not be described in detail, is the horizontal gradient mapping function, is the ray elevation angle, is the ray azimuth angle, The meanings of the formulas are the same, The meanings of the formulas are the same.

[0075] The calculation of the ray SWV value is completed, and the SWV matrix is constructed.

[0076] Step 4: Constructing the intercept matrix and the effective ray library according to the Beidou observation data at the tomographic time point. First, calculate the intercept of each ray in the voxels divided in the tomographic area, and then construct the intercept matrix of the ray. First, convert the geometric relationship between the receiver and the satellite to a three-dimensional voxel grid in the local rectangular coordinate system, express each signal path through the ray parameter equation, and calculate the intersection points of each voxel boundary. According to the start and end point parameters of the ray in the single voxel, the length of the straight line segment in the voxel is calculated, which is the intercept corresponding to the voxel. Then, the intercept matrix is constructed. Specifically, the present embodiment includes the following sub-steps:

[0077] (1) Establish local coordinate system and voxel grid: Convert the longitude, latitude and altitude of the Beidou receiver into a local rectangular coordinate system (obtained in step 1) and define a regular three-dimensional voxel grid in the coordinate system, and give the grid division points in each direction.

[0078] (2) Ray parameterization: Use the Beidou receiver position and the azimuth and elevation angle of the ray to construct the starting vector and the direction vector for each ray, and express it using the parameter equation:

[0079]

[0080] wherein is the ray parameter and is the part from the receiver to the satellite, is the ray direction vector, is the ray starting point, is the parameter equation of the ray.

[0081] (3) Pre-computation of ray and grid plane intersection points: For each ray, find the intersection with all x, y and z direction grid planes to obtain the corresponding parameter values , and eliminate parallel or numerical anomaly cases. This step obtains a set of candidate intersection point parameters.

[0082] (4) Determine valid intersection points voxel by voxel: For each voxel, select the corresponding 6 boundary planes, find the intersection point coordinates from the candidate parameters, and determine whether the intersection point falls on the surface of the voxel under the condition of numerical tolerance. According to whether the receiver is located inside the voxel, constrain the number of allowed intersection points (0, 1, 2) of the ray in the voxel, so that the actual number of intersection points satisfies the geometric topological relationship (not passing through, inside, passing through).

[0083] (5) Calculate the path length in the voxel: When the number of valid intersection points in a voxel is determined, the intercept length in the voxel can be calculated.

[0084] If there are two intersection point parameters and , then the intercept length in the voxel is:

[0085]

[0086] wherein is the intercept length in the voxel , is the ray direction vector, is the modulus of the ray direction vector, and are the parameters of the ray at the intersection point.

[0087] If the receiver is inside the voxel and there is only one intersection point then the length is:

[0088]

[0089] where, is the in-voxel intercept length , is the ray direction vector, is the magnitude of the ray direction vector, is the parameter of the ray at the intersection point.

[0090] If there is no intersection then:

[0091]

[0092] where, is the in-voxel intercept length .

[0093] This process is repeated for all rays and voxels to obtain the final intercept matrix.

[0094] After the parameter calculation of each ray is completed, the step of dividing the effective rays into multiple ordered balanced subsets includes: mixing and uniformly distributing all effective rays according to at least one attribute of the elevation angle, the azimuth angle, the receiver station to which the ray belongs, and the ray path length or the number of traversed voxels, so that each subset contains rays of different attribute ranges at the same time, and each subset is representative of the statistical characteristics of the entire ray set.

[0095] Specifically, due to the large error of low-elevation rays and the fact that rays emitted from the side of the tomography region contain water vapor information not belonging to the tomography region, this embodiment selects rays with a suitable elevation angle (greater than 15°) that emerge from the top or near the top of the tomography region (where the projection of the ray length in the tomography region in the zenith direction is greater than 9 km). Based on the geographical location of the receiving station, the ray elevation angle distribution, the sampling window, and other information, all effective rays are sequentially divided into multiple ordered equilibrium subsets. In this case, a sampling window of 10 minutes and a sampling rate of 2 minutes were selected. The sample was divided into 6 ordered subsets. Rays were evenly distributed among the subsets based on key geometric and observational attributes, ensuring each subset represented the whole as much as possible. First, bins were divided by elevation angle, then sectors by azimuth angle. Rays within each bin and sector were distributed alternately to ensure each subset contained rays from both high and low elevation angles and different directions. Simultaneously, a balanced distribution was achieved by station to prevent any subset from being overly concentrated on a few stations. The distribution of path length or the number of voxels traversed was also considered to prevent any subset from consisting entirely of short or long rays, thereby improving the stability and convergence consistency of iterative updates. Each subset contained most of the features of the whole, maintaining subset balance, and the intercept matrix was processed accordingly. The effective ray library composed of subsets, including parameters such as the SWV of the rays (obtained in step 3) and the intercept matrix, were stored in text format on the local computer. This completed the construction of the effective ray library composed of subsets and the establishment of the corresponding intercept matrix.

[0096] C. A fusion algebraic reconstruction algorithm is adopted, combined with dual-source relaxation factor constraints based on height prior information and iterative residual information, and iterative updates are performed using the effective ray library. The BeiDou tropospheric tomography matrix equations constructed from the iterative initial field, intercept matrix, and oblique path water vapor values ​​are solved, outputting the three-dimensional water vapor density inversion results, such as... Figure 3 , Figure 4 As shown, steps 5-9 are included:

[0097] Step 5: Predict the water vapor density at the tomography time point using the water vapor prediction model as the initial value for iteration. Calculate the recent water vapor density using radiosonde data from the radiosonde station, selecting a time period approximately one week before the tomography time point, with an interval of about 6-24 hours between each time step, for a total of 14 time steps of water vapor density data. Also use the PWV (Potential Water Vapor Value) of the HKSC station obtained through BeiDou data processing. Since the HKSC station is closest to the 45004 radiosonde station and located in the same tomography grid, HKSC station data is chosen for PWV calculation. The formula for PWV calculation using BeiDou data is as follows:

[0098]

[0099] In the formula, Atmospheric precipitable water , is the conversion factor, is the zenith wet delay , has the same meaning as in formula , has the same meaning as in formula (Step 3), has the same meaning as in formula .

[0100] The water vapor density of the first 14 time steps before the tomography time point, the annual accumulated day sin-cos mapping value at the tomography time point, and the PWV at the tomography time point of the HKSC site are input into the water vapor density prediction model based on LSTM fusion PWV. The model predicts the output of the water vapor density at the tomography time point, and then the predicted water vapor density is input into the tomography equation constructed in the next step as the initial value of iteration for iterative inversion.

[0101] Step 6: Constructing the Beidou troposphere tomography matrix equation. The intercept matrix of the ray, the SWV matrix, and the water vapor density matrix are numbered. According to the voxel division strategy (proposed in step 1), for the ray intercept matrix, each row represents all the intercept data of a ray, and each column represents a voxel grid. The grid that is not crossed by the ray has an intercept value of 0. For the SWV matrix, each row represents the SWV value of a ray. For the water vapor density matrix, all voxels are numbered by height and latitude expansion to form a row, and each column of the matrix represents a corresponding grid to each column of the intercept matrix. The initial value of the water vapor density matrix is calculated by the water vapor prediction model (obtained in step 5), and the water vapor density initial value of the same height layer is the same. The constructed tomography equation is as follows:

[0102]

[0103] In the formula, is the intercept matrix, is the water vapor density matrix, is the ray SWV value matrix.

[0104] Plus appropriate horizontal Gaussian empirical constraints and vertical exponential empirical constraints, which are common knowledge in the field and will not be described in detail, and the formula is as follows:

[0105]

[0106]

[0107] In the formula, is the horizontal empirical constraint matrix, is the vertical empirical constraint matrix, is the water vapor density matrix, has the same meaning as in formula The meanings are the same.

[0108] The main role of the experience constraint is to make certain adjustments to the grid water vapor density that is not penetrated by the ray, so that the water vapor density in the voxel that is not penetrated by the ray is closer to the true value.

[0109] Step 7: Calculate the high experience weight factor of the relaxation factor. After the tomography equation is constructed, the inversion is solved by a new optimization iteration method, and the experience weight constraint is assigned to the iterative relaxation factor of the high layer voxel according to the high experience water vapor model. Since the water vapor density of the high layer voxel is generally small and has universality and regularity relative to the water vapor distribution of the middle and low layers, an exponential normalization assignment model is selected. The specific assignment calculation formula is as follows:

[0110]

[0111] In the formula, is the high experience weight factor of the first voxel, is the altitude of the first voxel , is the altitude threshold , the relaxation factor weight of the voxel less than the altitude threshold is 1, is the altitude difference , that is, and absolute difference, is the water vapor elevation, which can be obtained by an empirical formula, which is common sense in the field. In this case, the altitude threshold is selected as 6km. After the initialization of the high experience weight factor of the relaxation factor is completed, the tomography equation is fused and iteratively reconstructed.

[0112] Step 8: Iteratively solve the tomography equation by the fusion algebraic reconstruction method of the double-source constrained relaxation factor. The initial value of the water vapor iteration is calculated by the water vapor prediction model (obtained in step 5).

[0113] As shown in Figure 4 , the fusion algebraic reconstruction method is specifically divided into two parts, one is the MART part, and the other is the SIRT part. MART is to iterate on each grid, and SIRT is to iterate on the whole. At the same time, MART iteration naturally maintains the non-negative property, which is suitable for water vapor density inversion. One MART iteration on each subset is called one iteration, and one SIRT iteration on the whole is performed after each Y T MART iteration, and SIRT is selected every 6 MART after traversing all subsets once to prevent local iteration problems. The iteration formulas of MART and SIRT are as follows:

[0114]

[0115]

[0116] wherein, is the water vapor density in voxel , is the voxel number, is the iteration number, is the water vapor density in the th voxel in the th iteration , is the water vapor density in the th voxel in the th iteration , is the ray number, is the ray number, is the voxel parameter number, is the intercept of the th ray in the th voxel in the intercept matrix , is the relaxation factor of the th voxel, is the SWV value of the th ray .

[0117] At the same time, the residual matrix of the overall ray SWV matrix and the inversion SWV matrix is calculated after each iteration, and the size of the iteration relaxation factor is controlled by the residual matrix, and the relaxation factor is updated after each iteration, and the high-level voxel also needs to meet the experience weight constraint. The representative parameters of the residual matrix need to be scaled or mapped, and here the average absolute residual calculated by the residual matrix is input into the function, and the value of the function is used as the result adjustment relaxation factor.

[0118] The final iteration relaxation factor and the calculation formula of the inversion reconstructed SWV value are:

[0119]

[0120]

[0121]

[0122]

[0123] wherein, is the relaxation factor of the th voxel, is the relaxation factor of the The height experience weight factor of the voxel, The representative parameter of the residual matrix, The weight factor obtained by mapping through the hyperbolic tangent function, The representative parameter of the residual matrix, the mean absolute residual value calculated from the residual matrix is selected, The number of rays, The ray number, The intercept of the first iteration of the first ray, The SWV value reconstructed by inversion of the first iteration of the first ray, , The number of voxel parameters, The voxel number, The intercept of the first ray in the first voxel in the intercept matrix, , The water vapor density in the first iteration of the first voxel, . The same meaning as in formula , , , , , , , The same meaning as in formula , .

[0124] The termination condition of iteration is selected when the RMSE change rate of the ray SWV matrix and the inversion SWV matrix is less than 1% or the iteration number is greater than a certain degree (in this case, it is selected to be greater than 500 rounds) The iteration is terminated.

[0125] Finally, combined with appropriate horizontal and vertical experience constraints, the water vapor density of the grid not passed by the ray is adjusted. Thus, the final three-dimensional water vapor density in the tomographic region is obtained.

[0126] Step 9: Save and output the inversion results.

[0127] After the above-mentioned Beidou tropospheric water vapor tomography algorithm processing, the problems that the existing method cannot provide water vapor density with high accuracy, strong reliability, good time and space resolution as the initial value of tomography iteration, resulting in poor tomography results and weak sensing of water vapor vertical movement can be effectively solved; the problems of large calculation amount and many iteration times caused by traversing all rays once can be effectively solved; the problems of insufficient tomography precision caused by excessive reverse update, single algorithm defects, inability to adapt to the degree of ill-conditioning of different grids, and poor iteration stability of the iterative inversion algorithm can be effectively solved. The accuracy, real-time performance and robustness of the Beidou tropospheric tomography are effectively improved, and the challenges in tomography can be better addressed, thereby providing more stable and efficient technical support for the tropospheric water vapor tomography.

[0128] The basic principles, main features and advantages of the present application are shown and described above. It should be understood by those skilled in the art that the present application is not limited by the above-mentioned embodiments, and the above-mentioned embodiments and descriptions in the specification are only preferred examples of the present application and are not intended to limit the present application. Without departing from the spirit and scope of the present application, various changes and improvements can be made to the present application, and these changes and improvements all fall within the scope of the claimed present application. The scope of protection of the present application is defined by the appended claims and their equivalents.

Claims

1. A BeiDou tropospheric water vapor tomography method based on LSTM and ray ordered subset constraint strategy, characterized in that, Includes the following steps: Step 1: Construct a water vapor density prediction model, and use the water vapor density prediction model to predict the water vapor density distribution at the target time of the chromatography as the initial field for iteration; The steps for constructing the water vapor density prediction model include: Historical radiosonde water vapor data were acquired, and water vapor density and atmospheric precipitable water volume (PWV) at multiple altitudes at multiple historical moments were calculated based on the radiosonde water vapor data. Construct a neural network model that includes an input feature layer, an LSTM sequence modeling layer, a ReLU activation function layer, a fully connected layer, and a loss regression layer; The neural network model is trained by taking a combination of historical water vapor density sequences, atmospheric precipitable water volume (PWV), and time-coded features as input, and the water vapor density at the target time as output, to obtain the water vapor density prediction model. Step 2: Based on BeiDou observation data, obtain the effective ray set and its corresponding oblique path water vapor value at the time of the tomographic target. Calculate the intercept of each ray in the preset three-dimensional voxel grid to construct the intercept matrix. Divide all effective rays into multiple ordered balanced subsets to construct an effective ray library containing the oblique path water vapor value, intercept matrix and ordered balanced subset of effective rays, and save it in text form. Step 3: Employ a fusion algebraic reconstruction algorithm, combining dual-source relaxation factor constraints based on altitude prior information and iterative residual information, and using the effective ray library for iterative updates, to solve the BeiDou tropospheric tomography matrix equation constructed with the iterative initial field, intercept matrix, and oblique path water vapor values, and output the three-dimensional water vapor density inversion result. The fusion algebraic reconstruction algorithm includes a periodic alternating iterative algorithm of the multiplicative algebraic reconstruction method MART and the joint iterative reconstruction algorithm SIRT.

2. The BeiDou tropospheric water vapor tomography method based on LSTM and ray ordered subset constraint strategy according to claim 1, characterized in that, The time encoding feature is the feature value obtained by performing sine and cosine function transformations on the yearly cumulative day of the target time.

3. The BeiDou tropospheric water vapor tomography method based on LSTM and ray ordered subset constraint strategy according to claim 1, characterized in that, The method for using the water vapor density distribution at the predicted target time of the chromatography as the initial field for iteration is as follows: N times the time point before the chromatography time point... W The water vapor density sequence at each time step, the atmospheric precipitable water volume (PWV) at each chromatographic time point, and the coded features at each chromatographic time point are combined and input into the water vapor density prediction model for prediction. The predicted water vapor density sequence is then used as the initial field for iteration.

4. The BeiDou tropospheric water vapor tomography method based on LSTM and ray ordered subset constraint strategy according to claim 1, characterized in that, Step 2, which involves obtaining the effective ray set and the water vapor value along the oblique path, includes: The total zenith delay (ZTD) of each receiver station is calculated based on BeiDou observation data, and the zenith static delay (ZHD) and zenith wet delay (ZWD) are calculated in combination with the meteorological parameters of the stations. The zenith wet delay (ZWD) is converted into atmospheric precipitable water vapor (PWV) using a conversion factor. Based on the mapping function and horizontal gradient model, and combined with the elevation and azimuth angles of each ray, the zenith wet delay (ZWD) is projected onto each oblique path to obtain the oblique path water vapor value (SWV) of each ray, and then the SWV matrix is ​​constructed.

5. The BeiDou tropospheric water vapor tomography method based on LSTM and ray ordered subset constraint strategy according to claim 1, characterized in that, The intercept matrix is ​​constructed as follows: First, the geometric relationship between the receiver and the satellite is transformed into a three-dimensional voxel mesh in a local rectangular coordinate system. Each signal path is represented by the ray parametric equation, and its intersection with the boundaries of each voxel is calculated. Based on the start and end parameters of the ray's traversal segment within a single voxel, the length of its straight line segment within that voxel is calculated, and this length is the intercept corresponding to that voxel. Then, the intercept matrix is ​​constructed.

6. The BeiDou tropospheric water vapor tomography method based on LSTM and ray ordered subset constraint strategy according to claim 1, characterized in that, The steps of dividing the effective ray into multiple ordered balanced subsets include: Based on at least one of the following attributes of the ray: elevation angle, azimuth angle, receiver station, and path length or number of voxels traversed, all effective rays are mixed and uniformly distributed so that each subset contains rays with different attribute ranges, and each subset is representative of the statistical characteristics of the overall ray set.

7. The BeiDou tropospheric water vapor tomography method based on LSTM and ray ordered subset constraint strategy according to claim 1, characterized in that, The periodic alternating iteration method is as follows: after completing a preset number of MART iterations on all ordered balanced subsets, a SIRT iteration based on all rays is performed.

8. The BeiDou tropospheric water vapor tomography method based on LSTM and ray ordered subset constraint strategy according to claim 1, characterized in that, The relaxation factor constraint based on high prior information is: Different relaxation factor weights are assigned to voxels at different altitudes. For voxels with an altitude higher than a preset threshold, their weights are reduced exponentially based on the difference in altitude between the voxels and the preset threshold. For voxels with an altitude lower than or equal to the preset threshold, their weights are the baseline values.

9. The BeiDou tropospheric water vapor tomography method based on LSTM and ray ordered subset constraint strategy according to claim 1, characterized in that, The relaxation factor constraint based on iterative residual information is: After each iteration, the residual between the currently inverted oblique path water vapor value and the actual oblique path water vapor value is calculated, and the relaxation factor for the next iteration is dynamically adjusted based on the statistics of the residual.

10. The BeiDou tropospheric water vapor tomography method based on LSTM and ray ordered subset constraint strategy according to claim 9, characterized in that, The relaxation factor is dynamically adjusted based on the statistics of the residuals as follows: the mean absolute residual value of the residuals is input into the hyperbolic tangent function for mapping, and the mapping result is used as the adjustment coefficient to adjust the relaxation factor.

11. The BeiDou tropospheric water vapor tomography method based on LSTM and ray ordered subset constraint strategy according to any one of claims 1 to 10, characterized in that, After iteratively solving the tomography matrix equation, empirical constraints in the horizontal and vertical directions are introduced to constrain the water vapor density values ​​in voxels that are not penetrated by rays.

12. The BeiDou tropospheric water vapor tomography method based on LSTM and ray ordered subset constraint strategy according to any one of claims 1 to 10, characterized in that, Rays with an elevation angle greater than a set angle threshold and whose path projection height in the zenith direction is greater than a set height threshold are selected as the effective rays.

Images