A method for extracting ionospheric scintillation factor based on DORIS dual-frequency observation

Abstract

This invention discloses a method for extracting ionospheric scintillation factors based on DORIS dual-frequency observations, belonging to the field of space weather monitoring and satellite navigation technology. By constructing a cycle slip detection and low elevation angle control method suitable for the DORIS system, the RODTI (Rate of Change in Total Electron Content) factor of the ionosphere is extracted. A two-dimensional grid background field of ROTI is constructed using GPS observation data from IGS stations, achieving spatiotemporal matching between DORIS and GPS scintillation factors. Furthermore, a support vector regression method considering latitudinal constraints is introduced to normalize and correct the RODTI amplitude, solving the problem of inconsistent scintillation factor amplitudes between different observation systems. This yields a more accurate RODTI with amplitudes consistent with ROTI, thus forming a complete method for extracting ionospheric scintillation factors based on DORIS dual-frequency observations. This method can reliably reflect the spatiotemporal distribution characteristics of ionospheric scintillation on a global scale, providing a new data source for constructing a multi-source ionospheric scintillation dense monitoring network covering oceans, polar regions, and other areas, and has broad applicability.

Description

TECHNICAL FIELD

[0001] The application is suitable for an ionospheric scintillation factor extraction method based on DORIS dual-frequency observation, and belongs to the technical field of space weather monitoring and satellite navigation. TECHNICAL BACKGROUND With the arrival of the 25th solar cycle peak year, ionospheric disturbances have significantly increased, leading to frequent ionospheric scintillation events. Ionospheric scintillation has become a key factor affecting the stability of global navigation satellite system (GNSS) positioning services. Global monitoring of ionospheric scintillation is an important means to eliminate its interference with GNSS. Therefore, global monitoring of ionospheric scintillation must be carried out to eliminate its impact. Generally, ionospheric scintillation monitoring requires the use of professional equipment, ionospheric scintillation monitoring receivers (ISMR), which can directly provide high-precision scintillation factors to quantitatively characterize the impact of scintillation on signals. However, ISMR equipment is expensive, and the number of stations worldwide is extremely sparse, making it difficult to support the construction of a dense monitoring network required for high-resolution global ionospheric scintillation models.

[0002] In view of the limitations of ISMR station distribution, current ionospheric scintillation monitoring mainly utilizes the widely distributed ground-based GNSS receiver network. Many researchers have inverted various scintillation factors based on GNSS observation data, and a relatively comprehensive scintillation monitoring has been formed globally. However, this ground-based monitoring mode has the defect of uneven geographical distribution of GNSS receiving stations. In land areas such as North America, Europe, and East Asia, the station network is dense, and relatively fine monitoring results can be provided; but in vast oceans, harsh polar regions, and sparsely populated areas, the stations are sparse or even missing. This uneven distribution leads to significant blind spots in global ionospheric scintillation monitoring, making it impossible to construct a global coverage and high-precision ionospheric scintillation monitoring network, and hindering the in-depth understanding of the physical mechanisms of scintillation events in key regions such as the South Atlantic Anomaly (SAA).

[0003] To break through the monitoring limitations of GNSS network in the ocean and other areas, new and complementary observation data sources are urgently needed. On this basis, DORIS is considered as an important supplementary observation method. In recent years, preliminary progress has been made in ionospheric research using DORIS data. Scholars have constructed a global vertical total electron content (VTEC) model using DORIS data through assimilation technology, and pointed out that it has a significant contribution to improving the modeling accuracy in the equatorial ionospheric anomaly (EIA) region. At the same time, relevant scholars have proposed a differential slant TEC analysis method by taking advantage of the characteristics that the frequency difference of DORIS signal is nearly 5 times that of GPS, and verified the feasibility and consistency of DORIS data in ionospheric model quality evaluation. However, existing researches almost focus on the extraction of large-scale parameters such as TEC using DORIS, and there is still a lack of systematic research on the description of rapid disturbance characteristics caused by small-scale irregularities and the inversion of ionospheric scintillation factors. SUMMARY

[0004] In view of the shortcomings of the existing research, a DORIS system suitable for cycle slip detection method is provided, the elevation control and cycle slip detection are used to improve the quality of the observation value, and then the ionospheric scintillation factor RODTI is constructed based on the DORIS dual-frequency carrier phase observation value. The Kriging interpolation method is used to generate the GPS ROTI two-dimensional grid background field, the spatio-temporal matching data of RODTI and GPS ROTI is obtained, and the support vector regression method considering the latitude constraint is used to correct the amplitude of RODTI, so that the RODTI with consistent amplitude and more accurate ROTI is obtained.

[0005] To achieve the above technical purposes, the application discloses an ionospheric scintillation factor extraction method based on DORIS dual-frequency observation, and the steps are as follows: a. The DORIS carrier phase observation sequence is subjected to data preprocessing of cut-off height angle rejection and cycle slip detection, and high-quality carrier phase data is obtained for ionospheric scintillation factor calculation; b. The ionospheric piercing point IPP is calculated in combination with the ground beacon station coordinates and the DORIS precise orbit file; c. The total electron content rate index RODTI is calculated based on the DORIS carrier phase observation value after data preprocessing; d. The ROTI two-dimensional grid background field is constructed and matched with RODTI in space and time; e. The amplitude of RODTI is corrected based on the ROTI reference; f. The consistency of RODTI and ROTI is analyzed based on the amplitude correction result of RODTI, and the global ionospheric scintillation monitoring capability of RODTI is verified.

[0006] Furthermore, combining the DORIS precise orbit files provided by the international DORIS service organization IDS with the DORIS dual-frequency carrier phase observation data RINEX, the DORIS carrier phase observation sequence underwent data preprocessing including cutoff elevation angle removal and cycle slip detection. The elevation angle magnitude for each epoch was calculated, and cycle slips were detected at each epoch. Epochs with cycle slips were removed to avoid introducing non-flickering phase abrupt changes during differential processing. The process is as follows: The elevation angle calculation process for each epoch is as follows: The position of the satellite equipped with the DORIS receiver is obtained using a precise orbit file; combined with the coordinates of the ground beacon station, the coordinate difference vector between the satellite and the ground station is obtained. : , In the formula Indicates the satellite's location. Indicates the coordinates of the ground beacon station; Coordinate difference vector The local northeast celestial coordinate components of the station center are obtained through the coordinate transformation matrix: , ; In the formula, N, E, and U are the north, east, and zenith components, respectively; T is the coordinate transformation matrix; and B and L are the latitude and longitude of the ground station, respectively. Based on this, the elevation angle of the satellite relative to the ground beacon station at each epoch can be calculated using the following formula. : , In the formula, Ele is the satellite elevation angle value. The cutoff elevation angle is set to 30°, and observations below 30° are discarded. The cycle slip detection process is as follows: The geometrically independent combination value for each epoch is calculated using the following formula. : , In the formula, i represents time. and These are dual-frequency carrier phase observations. and They are respectively and For the wavelength corresponding to the carrier frequency, the first-order difference value of the geometrically independent combination is calculated for each DORIS continuous observation arc segment using the following formula. The continuous observation arc refers to the effective observation sequence in the preprocessed DORIS dual-frequency observation data that is continuous in adjacent epochs, has no missing data, and does not cross cycle jump points. , Using third-order polynomials to evaluate first-order difference values Perform fitting to eliminate Trend term in the fit: , , In the formula The first-order difference value of the geometrically independent combination after fitting the third-order polynomial is given. ,…, These are the polynomial coefficients, which are determined by the least squares method. This is the starting epoch of the current continuous arc segment. This is the end epoch of the current continuous arc segment. These represent the observation times at each epoch of the current arc segment; Using the formula: Calculate the fitting residuals for each epoch. According to 3 Principle, when satisfied or When determining whether a cycle jump exists in the current epoch, among which... Let be the mean of the fitting residuals within the arc segment. denoted as the standard deviation of the fitting residuals within the arc segment.

[0007] Furthermore, the calculation process for the ionospheric penetration point (IPP) is as follows: Based on the DORIS observation data provided by the IDS center, using the ionospheric thin shell assumption model, and utilizing the GNSS station positions in the observation data, combined with the satellite positions in the precise ephemeris file, the IPP of the ionospheric penetration point (IPP) for all satellite line-of-sight directions of all GNSS stations is calculated according to the satellite-terrestrial geometry. , , , , , in, The zenith angle of the GNSS station. The elevation angle of the satellite relative to the GNSS station. The zenith angle of the puncture point. Let H be the Earth's average radius, taken as 6371 km, and let H be the ionospheric height, taken as 350 km. The latitude of the puncture point. For GNSS station latitude, This is the azimuth angle of the satellite relative to the GNSS station. Longitude of the puncture point This refers to the longitude of the GNSS station.

[0008] Furthermore, the calculation process for the total electron content change rate index RODTI is as follows: Common errors such as geometric distance, satellite and receiver clock bias are eliminated through a geometrically independent linear combination to obtain the pure ionospheric delay observation. , In the formula yes Geometrically independent linear combination equations of DORIS dual-frequency observations at time 10 ... and Here are the carrier phase observations, where λ1 and λ2 are the carrier frequencies, respectively. and The wavelength; To characterize the short-term variation of the total electron content (TEC) in the ionosphere, the rate of change (RODT) of TEC was first calculated based on a carrier phase combination sequence of adjacent epochs, with units of 1000 volts. The RODT between adjacent epochs is represented as: , In the formula t is the time interval between adjacent epochs, in minutes; assuming a sampling rate of 10s for extracting scintillation factors... ; Statistical analysis of the RODT sequences was performed using a 5-minute sliding window. The standard deviation of RODT within the window was defined as the DORIS ionospheric scintillation factor RODTI. , In the formula This indicates the epoch number corresponding to the current sliding window. Let N be the index of each epoch within the window, and let N be the number of epochs within the sliding window. Given a DORIS sampling interval of 10 seconds and a window length of 5 minutes, then N is 30. The average RODT value within time interval i. This represents the ionospheric scintillation factor value for each epoch.

[0009] Furthermore, spatiotemporal matching is performed between the constructed ROTI two-dimensional grid background field and RODTI, as follows: The ROTI (Relationship to Transient Regulator) of GNSS stations in each region was calculated using a sliding window statistical carrier phase change rate method. Under the assumption of a thin-shell ionospheric model, the effective height of the ionosphere was uniformly set to 350 km to determine the ionospheric penetration point (IPP) locations of each GNSS station at different epochs. Then, the irregular ROTI values ​​at the IPP points of each station were mapped onto a regular latitude and longitude grid using Kriging interpolation to generate a GPS ROTI two-dimensional grid background field covering the study area. The regular latitude and longitude grid is a set of grid nodes discretely generated within the study area according to the longitude and latitude intervals determined by the range of ionospheric penetration points in the GPS observation data, with a grid resolution set to 0.5°×0.5°. The generated ROTI two-dimensional grid is a two-dimensional continuous field composed of regular latitude and longitude grid points and their corresponding ROTI values ​​at each preset time level. Based on the constructed ROTI two-dimensional grid, the coordinates of the ionospheric puncture points and the observation time of each epoch on the continuous DORIS observation arc are extracted. The spatiotemporal bilinear interpolation method is used to obtain the ROTI values ​​corresponding to the DORIS observation times from the ROTI two-dimensional grid. The specific matching process is as follows: Let the IPP point of a certain epoch in DORIS be... In a ROTI 2D grid, the four nearest ROTI 2D grid points of a DORIS epoch IPP point are defined as follows: , , , The ROTI value corresponding to four adjacent grid points is , , , In a given time layer The IPP points are then calculated using spatial bilinear interpolation. ROTI: , In the formula , These represent the latitude and longitude of the IPP point corresponding to the DORIS observation epoch, respectively. , The grid spacing is between latitude and longitude. This represents the ROTI value obtained through spatial interpolation at a given time level; In the time dimension, for two adjacent time layers and Linear interpolation is performed to obtain ROTI values ​​that match RODTI after both spatial and temporal interpolation. : , In the formula and For adjacent time layers of the ROTI 2D grid, This indicates the actual time corresponding to the observed epoch; At each epoch of the DORIS trajectory arc, the corresponding ROTI value is extracted to obtain the matching data pair of RODTI and ROTI, forming the final matching dataset {ROTI, RODTI}, so as to compare and study the consistency of the two factors in monitoring ionospheric scintillation.

[0010] Furthermore, based on the matching dataset, the amplitude of RODTI is corrected using ROTI as a benchmark to reduce the systematic amplitude deviation between RODTI and ROTI caused by differences in observation frequency, signal propagation path, sampling method and spatial location, so that the corrected RODTI has better consistency with ROTI on a numerical scale. The RODTI amplitude correction method is as follows: a feature embedding machine learning model is constructed using support vector machine regression (SVR), a training sample set is built and input and output features are defined, the feature embedding machine learning model is trained using the training sample set to learn the mapping relationship between RODTI and ROTI, and finally the trained feature embedding machine learning model can be used to map and predict the input RODTI and the latitude to obtain the amplitude-corrected RODTI value. The constructed training sample set is represented as: , In the formula For the first The input feature vector of the matching sample is given by the first matching sample. DORIS scintillation factor corresponding to each matched sample and the The IOP latitude of each matched sample corresponding to the ionospheric puncture point. , For the first The regression target of each matched sample N is the total number of samples; Feature embedding machine learning model training parameters include regularization parameters Error tolerance bandwidth and nuclear scale ; The tradeoff between controlling model complexity and training error. Determine the range by which the model can ignore small errors. The smoothness of the kernel function and the model's sensitivity to input features are controlled; Bayesian optimization is used to automatically search for the optimal parameter combination to maximize the model's generalization ability; the hyperparameter search space is set as follows: , , During parameter search and model evaluation, 5-fold cross-validation is used, and the samples are randomly divided into five subsets for training and validation in turn, so as to make full use of the data and prevent overfitting. At the same time, the model stability is improved by 200 iterations of optimization while ensuring computational efficiency. Finally, the trained feature embedding machine learning model is used to map the input RODTI values ​​to the latitudes, resulting in normalized RODTI values: , In the formula This represents the regression mapping function between RODTI and ROTI obtained during training. This represents the normalized RODTI value.

[0011] Furthermore, the amplitude correction results of RODTI were verified globally to determine whether RODTI has the ability to monitor ionospheric scintillation. The verification method involved conducting spatial and temporal accuracy analyses at both the global and latitudinal levels. Simultaneously, the response capability of the extracted scintillation factor RODTI under geomagnetic storm events was evaluated. The process is as follows: a. Verification of spatial and temporal accuracy: We selected data from 250 GPS stations worldwide and observation data from 7 relatively active satellite platforms equipped with the DORIS system. We used amplitude-corrected RODTI and ROTI to conduct accuracy analysis on spatial, temporal, and typical geomagnetic storm events in global and latitudinal ranges. In spatial validation, matched samples formed over consecutive dates are statistically analyzed at both the global and latitudinal levels. The correlation and error characteristics of the amplitude-corrected RODTI and ROTI in different spatial regions are compared to verify the monitoring capability of the RODTI across global spatial distribution. The statistical indicators are as follows: R... 2 RODTI reflects the correlation between RODTI and RODTI; RMSE measures the overall dispersion of their residuals; Bias characterizes the magnitude and direction of systematic bias. 2 The higher the value, the smaller the RMSE, and the closer the bias is to 0, which proves that the monitoring results of the proposed RODTI are more reliable. In the time validation, to further verify the monitoring accuracy of the normalized RODTI on the time series, the daily average root mean square error RMSE and systematic bias Bias were calculated respectively. The calculated daily RMSE and Bias reflect the overall error magnitude and systematic bias of the normalized RODTI and ROTI on that day, respectively, and thus evaluate the stability and difference performance of the normalization method on the time series. b. To evaluate the response capability of RODTI under geomagnetic storm events, typical events during strong geomagnetic storms were selected as verification objects. The space weather background was determined by combining interplanetary magnetic field parameters and geomagnetic activity index. A time-continuous RODTI sequence was constructed by integrating observation data from multiple low-orbit satellites equipped with DORIS receivers. The RODTI retrieved from co-located GPS stations was used as a comparative reference. By analyzing the synchronous change characteristics of RODTI and RODTI during the main and recovery phases of the geomagnetic storm, the consistency and temporal reliability of RODTI's response to ionospheric disturbances caused by geomagnetic storm energy injection were evaluated, thereby verifying its monitoring capability under strong space weather conditions.

[0012] A computer device, characterized in that it includes a processor and a memory, the processor being electrically connected to the memory, the memory being used to store instructions and data, and the processor being used to execute an ionospheric scintillation factor extraction method based on DORIS dual-frequency observations.

[0013] Beneficial Effects: This invention proposes an ionospheric scintillation factor extraction method based on DORIS dual-frequency observations. By combining low elevation angle control with geometrically independent cycle slip detection, it effectively improves the availability and stability of carrier phase observation data, achieving reliable characterization of small-scale irregular volume disturbances in the ionosphere. This method fully utilizes the advantage of the large frequency interval in DORIS dual-frequency observations, enhancing the response capability to rapid ionospheric disturbance characteristics, expanding the application scope of the DORIS system in ionospheric scintillation monitoring, and providing a new data source for constructing a more refined global ionospheric scintillation monitoring network. Attached Figure Description

[0014] Figure 1 This is a flowchart illustrating the ionospheric scintillation factor extraction method based on DORIS dual-frequency observations according to the present invention. Detailed Implementation

[0015] The embodiments of the present invention will be further described below with reference to the accompanying drawings: like Figure 1 As shown, this invention discloses a method for extracting ionospheric scintillation factors based on DORIS dual-frequency observations, the steps of which are as follows: a. Combining the DORIS precise orbit files provided by IDS, and based on the DORIS dual-frequency carrier phase observation data provided by the IDS center, cutoff elevation angle removal and cycle slip detection processing are performed to improve the quality of the observation data.

[0016] Elevation angle calculation: First, the position of the satellite carrying the DORIS receiver is obtained using the precise orbit SP3 file. Then, combined with the coordinates of the ground beacon station, the coordinate difference vector between the satellite and the ground station is obtained. , In the formula Indicates the satellite's location. Indicates the coordinates of the ground beacon station. The coordinate difference vector pointing from the ground station to the satellite. Then, using a coordinate transformation matrix, the local northeast-sky coordinate component of the station center is obtained. , in, , In the formula, N, E, and U represent the north, east, and zenith direction components, respectively; T is the coordinate transformation matrix; and B and L are the latitude and longitude of the ground station, respectively. Based on this, the elevation angle of the satellite relative to the ground beacon station can be calculated. , In the formula, Ele represents the satellite elevation angle. The elevation angle value for each epoch is calculated using the above formula. The cutoff elevation angle is set to 30° to eliminate observations below 30°.

[0017] Cycle slip detection: First, calculate the geometry-free (GF) combination value for each epoch. , In the formula and These are dual-frequency carrier phase observations. This corresponds to the wavelength of the carrier frequency. Subsequently, calculations are performed for each consecutive observation arc segment. The first difference of the combination, , In the formula This is the first-order difference value of the GF combination. To eliminate... The trend term in the equation is fitted to the first-order difference value using a third-order polynomial. The formula for polynomial fitting is as follows. , , In the formula After polynomial fitting Combining first-order difference values, ,…, These are the polynomial coefficients, which are determined by the least squares method. This is the starting epoch of the current continuous arc segment. Let t be the end epoch of the current continuous arc segment, and t be the observation time of each epoch of the current arc segment.

[0018] Based on this, the fitting residuals for each epoch are calculated. , According to 3 In principle, when or The time period determines whether a cycle jump exists in the current epoch. and These represent the mean and standard deviation of the fitted residuals within each arc segment. After identifying cycle slips in each arc segment, the cycle slip repair method for DORIS phase observations is not discussed at this stage; instead, the epochs where cycle slips occurred are directly removed to avoid introducing non-flickering phase abrupt changes during differential processing.

[0019] b. Calculate the ionospheric puncture point (IPP) based on IDS dual-frequency carrier phase observations, and calculate the Rate of DORIS TEC Index (RODTI).

[0020] Calculation of Ionospheric Penetration Point (IPP): Based on DORIS observation data provided by the IDS Center, using the thin-shell ionospheric shell assumption model, and setting the ionospheric height to 350 km, the IPP of all stations along all satellite line-of-sight directions was calculated according to satellite-terrestrial geometry, utilizing the station locations in the observation data and combining them with satellite locations in the precise ephemeris files. , , , , , in, The zenith angle of the station. The satellite's elevation angle relative to the station. The zenith angle of the puncture point. Let H be the Earth's average radius, taken as 6371 km, and let H be the ionospheric height, taken as 350 km. The latitude of the puncture point. The latitude of the measuring station. This is the azimuth angle of the satellite relative to the station. Longitude of the puncture point The longitude of the station; c. Calculate the RODTI (Rate of Change in Total Electron Content): The RODTI, characterizing ionospheric scintillation, is calculated using carrier phase observations from DORIS data. First, common errors such as geometric distance, satellite and receiver clock bias are eliminated through a geometrically independent linear combination to obtain the pure ionospheric delay observation. , In the formula, L(i) is the geometrically independent linear combination equation of the DORIS dual-frequency observations at time i. and λ1 and λ2 are the carrier phase observations, and λ1 and λ2 are the wavelengths of the corresponding carrier frequencies.

[0021] The rate of change of TEC (RODT) is calculated using time difference, with units of 1000 volts per kilometre var. tekt (TDT). It can reflect the instantaneous intensity of ionospheric disturbances: , In the formula t is the time interval between adjacent epochs (unit: min). The scintillation factor is extracted at a sampling rate of 10s. t = 10 / 60 ≈ 0.1667 min.

[0022] To capture the perturbation characteristics caused by ionospheric irregularities in the RODT sequence, a 5-minute sliding window was used to calculate its standard deviation, and the scintillation factor was defined as RODTI. , In the formula, N is the number of epochs within the sliding window. is the average RODT value within time interval i.

[0023] d. The accuracy of RODTI is verified and analyzed based on the Rate of TECIndex (ROTI) calculated from global GPS observation data with a sampling interval of 30 seconds provided by IGS. First, a two-dimensional grid background field for ROTI is constructed using GPS observation data from IGS stations, and then spatiotemporally matched with RODTI. The method for calculating the ionospheric puncture point (IPP) and ROTI using GPS is the same as that used in DORIS.

[0024] Constructing a GPS ROTI 2D grid background field: First, the ROTI of GNSS stations in each region was calculated using the same method. Under the assumption of a thin-shell ionospheric model, the effective ionospheric height was uniformly set to 350 km to determine the ionospheric penetration point (IPP) locations of each station at different epochs. Subsequently, the irregular ROTI values ​​at the IPPs of each station were mapped to a regular grid using Kriging interpolation, thereby generating a ROTI 2D grid covering the entire region.

[0025] To achieve spatiotemporal matching between RODTI and ROTI: The coordinates of the ionospheric puncture point (IPP) and the observation time of each epoch on the continuous observation arc of DORIS are extracted. A spatiotemporal bilinear interpolation method is used to obtain the ROTI values ​​corresponding to the DORIS observation times from the grid. Let the IPP position of a certain epoch in DORIS be... In a ROTI grid, the four neighboring grid points are defined as follows: , , , The corresponding ROTI value is , , , Then, in a given time layer... The ROTI of the IPP point can be calculated using spatial bilinear interpolation. , In the formula , These represent the latitude and longitude of the IPP corresponding to the DORIS observation epoch, respectively. , The grid spacing is between latitude and longitude. This represents the ROTI value obtained through spatial interpolation at a given time level.

[0026] Furthermore, in the time dimension, for two adjacent time layers and Linear interpolation yields the ROTI value at the actual time step. , In the formula and For adjacent time layers of the ROTI grid, This represents the actual time corresponding to the observed epoch. The ROTI value is obtained by spatial and temporal interpolation and matches the RODTI.

[0027] By using the above spatial and temporal dual interpolation, the corresponding ROTI value can be extracted at each epoch of the DORIS trajectory arc, obtaining matching data pairs of RODTI and ROTI, forming a matching dataset {ROTI, RODTI}, so as to compare and study the consistency of the two factors in monitoring ionospheric scintillation.

[0028] e. Given the significant difference between the DORIS system frequency and the GPS system frequency, the amplitude range of RODTI is larger than that of ROTI. The amplitude of RODTI is corrected by using ROTI as a reference for further verification.

[0029] The RODTI amplitude correction method proposes a machine learning-based RODTI normalization method that uses ROTI as a reference and incorporates latitudinal constraints. This method embeds latitude as a key input feature into the machine learning model, guiding it to learn the characteristics of ionospheric perturbation variation with latitude, thereby enhancing the spatial dependence of the normalization result and helping to avoid unreasonable systemic biases. This method selects Support Vector Regression (SVR) as the machine learning tool to implement this normalization method. The specific implementation process is as follows: First, a training sample set is constructed and input / output features are defined. The RODTI of each observed arc segment is combined with its corresponding latitude information to form the input feature vector. The RODTI at the same time and location is used as the regression target, thus forming a supervised learning sample set. , In the formula For the input feature vector, The regression target is N, which is the total number of samples.

[0030] Then, the constructed sample set is input into the SVR model for training to learn the mapping relationship between RODTI and ROTI. Key parameters for SVR model training include the regularization parameter. Error tolerance bandwidth and nuclear scale .in The tradeoff between controlling model complexity and training error. Determine the range by which the model can ignore small errors. This paper controls the smoothness of the kernel function and the model's sensitivity to input features. Bayesian optimization is used to automatically search for the optimal parameter combination to maximize the model's generalization ability. The search space for hyperparameters is set as follows: , , During parameter search and model evaluation, 5-fold cross-validation was used, randomly dividing the samples into five subsets for training and validation in turn, in order to make full use of the data and prevent overfitting. At the same time, the model stability was improved by 200 iterations of optimization while ensuring computational efficiency.

[0031] Finally, using the trained SVR model, the input RODTI is mapped to the latitude to predict, and the normalized RODTI value is obtained: , In the formula This represents the regression mapping function between RODTI and ROTI obtained during training. This represents the normalized RODTI value.

[0032] f. Select data from 250 GPS stations worldwide and observation data from 7 relatively active satellite platforms equipped with the DORIS system. Use amplitude-corrected RODTI and ROTI to perform accuracy analysis on spatial, temporal, and typical geomagnetic storm events in global and latitudinal ranges.

[0033] In spatial validation, matched samples formed over consecutive dates are statistically analyzed at both the global and latitudinal levels. The correlation and error characteristics of the amplitude-corrected RODTI and ROTI in different spatial regions are compared to verify the monitoring capability of the RODTI across global spatial distribution. The statistical indicators are as follows: R... 2 The R² reflects the correlation between RODTI and ROTI. RMSE is used to measure the overall dispersion of the residuals of the two. Bias characterizes the magnitude and direction of systematic bias. If the R² is higher, the RMSE is smaller, and the Bias is closer to 0, it proves that the monitoring results of the proposed RODTI are more reliable. In the time-series validation, to further verify the monitoring accuracy of the normalized RODTI over time, the daily average root mean square error (RMSE) and systematic bias (Bias) were calculated. The calculated daily RMSE and Bias reflect the overall error magnitude and systematic bias of the normalized RODTI and ROTI on that day, respectively, thus evaluating the stability and variability performance of the normalization method over time. The evaluation method is as follows: First, a matching dataset of GPS and DORIS is constructed using the method proposed in claim 2, and then the coefficient of determination (R²) is used... 2 RODTI is evaluated using root mean square error (RMSE) and systematic bias (Bias). 2 The R² value reflects the correlation between RODTI and ROTI. RMSE measures the overall dispersion of their residuals, while Bias characterizes the magnitude and direction of systematic bias. A higher R², a smaller RMSE, and a Bias closer to 0 indicate more reliable monitoring results from the proposed RODTI.

[0034] Analysis of Typical Geomagnetic Storm Events: To evaluate the response capability of RODTI under geomagnetic storm events, typical events during strong geomagnetic storms were selected as verification objects. Space weather background was determined by combining interplanetary magnetic field parameters and geomagnetic activity indices. Based on this, a time-continuous RODTI sequence was constructed by integrating observation data from multiple low-Earth orbit satellites equipped with DORIS receivers, and the RODTI retrieved from co-located GPS stations was used as a comparative reference. By analyzing the synchronous variation characteristics of RODTI and ROTI during the main and recovery phases of the geomagnetic storm, the consistency and temporal reliability of RODTI's response to ionospheric disturbances caused by geomagnetic storm energy injection were evaluated, thereby verifying its monitoring capability under strong space weather conditions.

[0035] The above description is merely one embodiment of the present invention and is not intended to limit the present invention. Any minor modifications, equivalent substitutions, and improvements made to the above embodiment based on the technical essence of the present invention should be included within the protection scope of the present invention.

Claims

1. A method for extracting ionospheric scintillation factors based on DORIS dual-frequency observations, characterized in that, The steps are as follows: a. Perform cutoff elevation angle removal and cycle slip detection on the DORIS carrier phase observation sequence to obtain high-quality carrier phase data for calculating the ionospheric scintillation factor; b. Calculate the Ionospheric Penetration Point (IPP) by combining the coordinates of ground beacon stations with the DORIS precision orbit file; c. Calculate the total electron content change rate index RODTI based on the DORIS carrier phase observations after data preprocessing; d. Construct a ROTI 2D grid background field and perform spatiotemporal matching with RODTI; e. Adjust the amplitude of RODTI based on ROTI; f. Analyze the consistency between RODTI and ROTI based on the amplitude correction results of RODTI, and verify the ability of RODTI to monitor global ionospheric scintillation.

2. The method for extracting ionospheric scintillation factors based on DORIS dual-frequency observations according to claim 1, characterized in that, Combining the DORIS precise orbit files provided by the international DORIS service organization IDS with the DORIS dual-frequency carrier phase observation data RINEX, the DORIS carrier phase observation sequence underwent data preprocessing including cutoff elevation angle removal and cycle slip detection. The elevation angle magnitude at each epoch was calculated, and cycle slips were detected at each epoch. Epochs with cycle slips were removed to avoid introducing non-flickering phase abrupt changes during differential processing. The process is as follows: The elevation angle calculation process for each epoch is as follows: The position of the satellite equipped with the DORIS receiver is obtained using a precise orbit file; combined with the coordinates of the ground beacon station, the coordinate difference vector between the satellite and the ground station is obtained. : ， In the formula Indicates the satellite's location. Indicates the coordinates of the ground beacon station; Coordinate difference vector The local northeast celestial coordinate components of the station center are obtained through the coordinate transformation matrix: ， ； In the formula, N, E, and U are the north, east, and zenith components, respectively; T is the coordinate transformation matrix; and B and L are the latitude and longitude of the ground station, respectively. Based on this, the elevation angle of the satellite relative to the ground beacon station at each epoch can be calculated using the following formula. : ， In the formula, Ele is the satellite elevation angle value. The cutoff elevation angle is set to 30°, and observations below 30° are discarded. The cycle slip detection process is as follows: The geometrically independent combination value for each epoch is calculated using the following formula. : ， In the formula, i represents time. and These are dual-frequency carrier phase observations. and They are respectively and For the wavelength corresponding to the carrier frequency, the first-order difference value of the geometrically independent combination is calculated for each DORIS continuous observation arc segment using the following formula. The continuous observation arc refers to the effective observation sequence in the preprocessed DORIS dual-frequency observation data that is continuous in adjacent epochs, has no missing data, and does not cross cycle jump points. ， Using third-order polynomials to evaluate first-order difference values Perform fitting to eliminate Trend term in the fit: ， ， In the formula The first-order difference value of the geometrically independent combination after fitting the third-order polynomial is given. ,…, These are the polynomial coefficients, which are determined by the least squares method. This is the starting epoch of the current continuous arc segment. This is the end epoch of the current continuous arc segment. These represent the observation times at each epoch of the current arc segment; Using the formula: Calculate the fitting residuals for each epoch. According to 3 Principle, when satisfied or When determining whether a cycle jump exists in the current epoch, among which... Let be the mean of the fitting residuals within the arc segment. denoted as the standard deviation of the fitting residuals within the arc segment.

3. The method for extracting ionospheric scintillation factors based on DORIS dual-frequency observations according to claim 2, characterized in that, The calculation process for the ionospheric penetration point (IPP) is as follows: Based on the DORIS observation data provided by the IDS center, using the thin ionospheric shell assumption model, and utilizing the GNSS station positions in the observation data, combined with the satellite positions in the precise ephemeris file, the IPP of the ionospheric penetration point (IPP) for all satellite line-of-sight directions of all GNSS stations is calculated according to the satellite-terrestrial geometry. ， ， ， ， ， in, The zenith angle of the GNSS station. The elevation angle of the satellite relative to the GNSS station. The zenith angle of the puncture point. Let H be the Earth's average radius, taken as 6371 km, and let H be the ionospheric height, taken as 350 km. The latitude of the puncture point. For GNSS station latitude, This is the azimuth angle of the satellite relative to the GNSS station. Longitude of the puncture point This refers to the longitude of the GNSS station.

4. The method for extracting ionospheric scintillation factors based on DORIS dual-frequency observations according to claim 3, characterized in that, The calculation process of the RODTI (Rate of Change in Total Electron Content) index is as follows: Common errors such as geometric distance, satellite and receiver clock bias are eliminated through a geometrically independent linear combination to obtain the pure ionospheric delay observation. ， In the formula yes Geometrically independent linear combination equations of DORIS dual-frequency observations at time 10 ... and Here are the carrier phase observations, where λ1 and λ2 are the carrier frequencies, respectively. and The wavelength; To characterize the short-term variation of the total electron content (TEC) in the ionosphere, the rate of change (RODT) of TEC was first calculated based on the carrier phase combination sequence of adjacent epochs, with units of 1000 m / s. The RODT between adjacent epochs is represented as: ， In the formula t is the time interval between adjacent epochs, in minutes; assuming a sampling rate of 10s for extracting scintillation factors... ; Statistical analysis of the RODT sequences was performed using a 5-minute sliding window. The standard deviation of RODT within the window was defined as the DORIS ionospheric scintillation factor RODTI. ， In the formula This indicates the epoch number corresponding to the current sliding window. Let N be the index of each epoch within the window, and let N be the number of epochs within the sliding window. Given a DORIS sampling interval of 10 seconds and a window length of 5 minutes, then N is 30. The average RODT value within time interval i. This represents the ionospheric scintillation factor value for each epoch.

5. The method for extracting ionospheric scintillation factors based on DORIS dual-frequency observations according to claim 4, characterized in that, Spatiotemporal matching was performed between the constructed ROTI 2D grid background field and RODTI, as follows: The ROTI (Relationship to Transient Regulator) of GNSS stations in each region was calculated using a sliding window statistical carrier phase change rate method. Under the assumption of a thin-shell ionospheric model, the effective height of the ionosphere was uniformly set to 350 km to determine the ionospheric penetration point (IPP) locations of each GNSS station at different epochs. Then, the irregular ROTI values ​​at the IPP points of each station were mapped onto a regular latitude and longitude grid using Kriging interpolation to generate a GPS ROTI two-dimensional grid background field covering the study area. The regular latitude and longitude grid is a set of grid nodes discretely generated within the study area according to the longitude and latitude intervals determined by the range of ionospheric penetration points in the GPS observation data, with a grid resolution set to 0.5°×0.5°. The generated ROTI two-dimensional grid is a two-dimensional continuous field composed of regular latitude and longitude grid points and their corresponding ROTI values ​​at each preset time level. Based on the constructed ROTI two-dimensional grid, the coordinates of the ionospheric puncture points and the observation time of each epoch on the continuous DORIS observation arc are extracted. The spatiotemporal bilinear interpolation method is used to obtain the ROTI values ​​corresponding to the DORIS observation times from the ROTI two-dimensional grid. The specific matching process is as follows: Let the IPP point of a certain epoch in DORIS be... In a ROTI 2D grid, the four nearest ROTI 2D grid points of a DORIS epoch IPP point are defined as follows: , , , The ROTI value corresponding to four adjacent grid points is , , , In a given time layer The IPP points are then calculated using spatial bilinear interpolation. ROTI: ， In the formula , These represent the latitude and longitude of the IPP point corresponding to the DORIS observation epoch, respectively. , The grid spacing is between latitude and longitude. This represents the ROTI value obtained through spatial interpolation at a given time level; In the time dimension, for two adjacent time layers and Linear interpolation is performed to obtain ROTI values ​​that match RODTI after both spatial and temporal interpolation. : ， In the formula and For adjacent time layers of the ROTI 2D grid, This indicates the actual time corresponding to the observed epoch; At each epoch of the DORIS trajectory arc, the corresponding ROTI value is extracted to obtain the matching data pair of RODTI and ROTI, forming the final matching dataset {ROTI, RODTI}, so as to compare and study the consistency of the two factors in monitoring ionospheric scintillation.

6. The method for extracting ionospheric scintillation factors based on DORIS dual-frequency observations according to claim 5, characterized in that, Based on the matching dataset, the amplitude of RODTI is corrected by using ROTI as a benchmark to reduce the systematic amplitude deviation between RODTI and ROTI caused by differences in observation frequency, signal propagation path, sampling method and spatial location, so that the corrected RODTI has better consistency with ROTI on a numerical scale. The RODTI amplitude correction method is as follows: a feature embedding machine learning model is constructed using support vector machine regression (SVR), a training sample set is built and input and output features are defined, the feature embedding machine learning model is trained using the training sample set to learn the mapping relationship between RODTI and ROTI, and finally the trained feature embedding machine learning model can be used to map and predict the input RODTI and the latitude to obtain the amplitude-corrected RODTI value. The constructed training sample set is represented as follows: ， In the formula For the first The input feature vector of the matching sample is derived from the first matching sample. DORIS scintillation factor corresponding to each matched sample and the Each matched sample corresponds to the IPP latitude of the ionospheric puncture point. , For the first The regression target of each matched sample N is the total number of samples; Feature embedding machine learning model training parameters include regularization parameters Error tolerance bandwidth and nuclear scale ; The tradeoff between controlling model complexity and training error. Determine the range by which the model can ignore small errors. The smoothness of the kernel function and the model's sensitivity to input features are controlled; Bayesian optimization is used to automatically search for the optimal parameter combination to maximize the model's generalization ability; the hyperparameter search space is set as follows: , , During parameter search and model evaluation, 5-fold cross-validation is used, and the samples are randomly divided into five subsets for training and validation in turn, so as to make full use of the data and prevent overfitting. At the same time, the model stability is improved by 200 iterations of optimization while ensuring computational efficiency. Finally, the trained feature embedding machine learning model is used to map the input RODTI values ​​to the latitudes, resulting in normalized RODTI values: ， In the formula This represents the regression mapping function between RODTI and ROTI obtained during training. This represents the normalized RODTI value.

7. The method for extracting ionospheric scintillation factors based on DORIS dual-frequency observations according to claim 6, characterized in that, The global accuracy of the RODTI amplitude correction results was verified to determine whether RODTI has the ability to monitor ionospheric scintillation. The verification method involved conducting spatial and temporal accuracy analyses at both global and latitudinal levels. Simultaneously, the response capability of the extracted scintillation factor RODTI under geomagnetic storm events was evaluated. The process is as follows: a. Verification of spatial and temporal accuracy: We selected data from 250 GPS stations worldwide and observation data from 7 relatively active satellite platforms equipped with the DORIS system. We used amplitude-corrected RODTI and ROTI to conduct accuracy analysis on spatial, temporal, and typical geomagnetic storm events in global and latitudinal ranges. In spatial validation, matched samples formed over consecutive dates are statistically analyzed at both the global and latitudinal levels. The correlation and error characteristics of the amplitude-corrected RODTI and ROTI in different spatial regions are compared to verify the monitoring capability of the RODTI across global spatial distribution. The statistical indicators are as follows: R... 2 RODTI reflects the correlation between RODTI and RODTI; RMSE measures the overall dispersion of their residuals; Bias characterizes the magnitude and direction of systematic bias. 2 The higher the value, the smaller the RMSE, and the closer the bias is to 0, which proves that the monitoring results of the proposed RODTI are more reliable. In the time validation, to further verify the monitoring accuracy of the normalized RODTI on the time series, the daily average root mean square error RMSE and systematic bias Bias were calculated respectively. The calculated daily RMSE and Bias reflect the overall error magnitude and systematic bias of the normalized RODTI and ROTI on that day, respectively, and thus evaluate the stability and difference performance of the normalization method on the time series. b. To evaluate the response capability of RODTI under geomagnetic storm events, typical events during strong geomagnetic storms were selected as verification objects. The space weather background was determined by combining interplanetary magnetic field parameters and geomagnetic activity index. A time-continuous RODTI sequence was constructed by integrating observation data from multiple low-orbit satellites equipped with DORIS receivers. The RODTI retrieved from co-located GPS stations was used as a comparative reference. By analyzing the synchronous change characteristics of RODTI and RODTI during the main and recovery phases of the geomagnetic storm, the consistency and temporal reliability of RODTI's response to ionospheric disturbances caused by geomagnetic storm energy injection were evaluated, thereby verifying its monitoring capability under strong space weather conditions.

8. A computer device, characterized in that, It includes a processor and a memory, the processor being electrically connected to the memory, the memory being used to store instructions and data, and the processor being used to execute the ionospheric scintillation factor extraction method based on DORIS dual-frequency observation as described in any one of claims 1-7.

Images