Volcanic eruption monitoring method by fusing data of Fengyun satellite and COSMIC ionospheric occultation

Abstract

The application discloses a volcanic eruption monitoring method fusing Fengyun satellite and COSMIC ionosphere occultation data, first, data is collected and preprocessed, the preprocessed data is fitted through spherical harmonics, and vertical total electron content of each grid point is obtained through inversion; abnormal state and normal state are defined according to solar radiation flux index under the grid, and dynamic threshold values of the abnormal state and the normal state are obtained; real-time data is collected, global grid vertical total electron content real-time values are obtained through inversion, are used for judging as suspected abnormal disturbance signals and recording relevant space-time information; vertical total electron content data of a preset area in a time window corresponding to the suspected abnormal disturbance signal is extracted, grid proportion is calculated through conditional calculation, and the abnormal disturbance signal is determined as solar or geomagnetic activity interference according to the proportion, otherwise, the abnormal disturbance signal is determined as a local volcanic eruption. The method can realize efficient, rapid, global coverage and high-precision volcanic eruption monitoring, and provides timely and reliable decision support for disaster prevention and reduction departments.

Description

Technical Field

[0001] This invention relates to the field of satellite remote sensing monitoring and disaster early warning technology, specifically a method for monitoring volcanic eruptions that integrates data from the Fengyun satellite and COSMIC ionospheric occultation. Background Technology

[0002] Volcanic eruptions are among the most devastating natural disasters affecting human life and the environment. They are often accompanied by seismic activity, causing not only crustal displacement but also impacting the ionosphere and upper atmosphere. Volcanic eruptions trigger a series of atmospheric pressure waves, which, through energy exchange with the ionosphere, disturb the upper atmosphere above the epicenter, further disrupting the ionosphere and affecting the normal operation of satellite navigation and communication systems. Statistics show that approximately 50-70 volcanoes erupt globally each year, with strong eruptions potentially causing global ionospheric disturbances that severely impact cross-regional communication and navigation. Therefore, timely and accurate volcanic eruption monitoring provides crucial support for disaster prevention and mitigation, reducing casualties and economic losses.

[0003] Currently, the main method for monitoring volcanic eruptions based on Total Electron Content (TEC) of the ionosphere is time-series analysis using ground-based GNSS data. This method calculates TEC using dual-frequency GNSS observations and detects ionospheric anomalies by analyzing changes in the TEC time series, thereby monitoring volcanic eruptions. Although this method has the advantage of high spatiotemporal resolution, the large volume of ground-based GNSS observation data and the complex multi-step data processing procedures result in a long overall processing time, which to some extent limits its effectiveness in large-scale monitoring and rapid response to volcanic eruptions. Therefore, researching efficient, rapid-response, wide-coverage, and high-spatiotemporal-resolution methods for volcanic eruption monitoring, providing strong support for natural disaster monitoring departments, is a pressing and challenging issue that needs to be addressed. Summary of the Invention

[0004] The purpose of this invention is to overcome the shortcomings of the prior art and propose a method for monitoring volcanic eruptions that integrates data from Fengyun satellites and COSMIC ionospheric occultation data. This method solves the problems of existing ground-based GNSS-based monitoring methods, such as complex and time-consuming data processing, difficulty in achieving rapid global response, and high false alarm rates due to interference from solar and geomagnetic activities. By integrating occultation data from COSMIC-1 and Fengyun-3C satellites, a 1°×1° gridded interpolation of the global ionosphere is performed using spherical harmonic functions. Combined with the geomagnetic index (Kp) and solar radiation flux (10.7 cm / 2800 MHz, F10.7) provided by the SEPC website as indicators of geomagnetic and solar activity intensity, a global ionospheric normal state background model and an ionospheric anomalous state background model are constructed, laying a benchmark for accurate extraction of disturbance signals. Based on occultation data from Fengyun-3D to Fengyun-3G satellites and COSMIC-2, VTEC is calculated, and the rate of change is calculated in real time. Combined with dynamic thresholds, this enables rapid identification of suspected anomalous signals, improving monitoring timeliness. The Kriging interpolation algorithm is used to model the VTEC of the spatial region with vertical total electron content above 60° north and south latitude. By calculating the rate of change, comparing it with thresholds, and analyzing regional disturbance indices, the interference of volcanic eruptions with solar and geomagnetic activity is accurately distinguished, reducing the false alarm rate. In summary, this invention aims to achieve efficient, rapid, globally covered, and highly accurate volcanic eruption monitoring, providing timely and reliable decision support for disaster prevention and mitigation departments.

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

[0006] This invention proposes a method for monitoring volcanic eruptions by fusing data from the Fengyun satellite and the COSMIC ionospheric occultation, comprising the following steps:

[0007] Step 1: Grid the globe and mark the grids containing active and dormant volcanoes. Collect historical ionospheric occultation data from COSMIC-1 and Fengyun-3C satellites during periods without VEI4 or higher volcanic eruptions. After preprocessing the data, remove outliers. Fit the preprocessed data with a spherical harmonic function to invert and obtain the total vertical electron content of each grid point.

[0008] Step 2: Define abnormal and normal states based on the solar radiation flux index under the grid, and obtain the dynamic thresholds for abnormal and normal states by calculating the standard deviation of the hourly change rate of the vertical total electron content at each grid point under the two states.

[0009] Step 3: Real-time acquisition of ionospheric occultation data from Fengyun-3D to Fengyun-3G satellites and COSMIC-2 satellite. After preprocessing and spherical harmonic function fitting, the real-time value of the vertical total electron content of the global grid is obtained by inversion. The data update frequency is 60 minutes / time. The rate of change of the vertical total electron content of each grid per hour is calculated. If the rate of change exceeds the dynamic threshold of the normal state and the grid is a marked volcanic area grid, it is judged as a suspected abnormal disturbance signal and the relevant spatiotemporal information is recorded.

[0010] Step 4: Extract the vertical total electron content data of the preset area within the time window corresponding to the suspected abnormal disturbance signal, and use the Kriging interpolation algorithm to construct the vertical total electron content spatial region of the grid within the preset area. Calculate the overall disturbance index and the proportion of grids exceeding the threshold within the time window based on the vertical total electron content spatial region. Determine whether it is solar or geomagnetic activity interference based on the proportion; otherwise, determine whether it is a local volcanic eruption.

[0011] Step 5: Add the collected and verified volcanic eruption event data and corresponding ionospheric observation data to the historical database and update the dynamic threshold.

[0012] Preferably, in step 1, during the gridding process, the globe is divided into 1°×1° grids and grids containing active volcanoes and dormant volcanoes are marked.

[0013] Preferably, in step 1, the preprocessing includes determining the quality of the electron profile based on the average deviation and the noise level factor and removing abnormal data, wherein the average deviation of the effective electron density profile is 0-1.5 and the noise level factor is less than 0.01.

[0014] Preferably, in step 1, the method for calculating the total vertical electron content is as follows: the electron density profile in the preprocessed ionospheric occultation data is interpolated using height splines, and the electron density in the 110-750km altitude range is numerically integrated to obtain the total oblique electron content. The total oblique electron content is then converted into a scatter plot of the total vertical electron content using a mapping function. The scatter plot of the total vertical electron content is then fitted with a spherical harmonic function to obtain the total vertical electron content in a 1°×1° grid.

[0015] Preferably, the mapping function is the cosine of the total oblique electron content multiplied by the zenith angle of the penetration point.

[0016] Preferably, in step 2, dates with a solar radiation flux index of not less than 150 sfu and an average geomagnetic index of not less than 4 are defined as ionospheric anomalous states, while the rest are normal states.

[0017] Preferably, in step 2, the rate of change of the total vertical electron content per hour is the total vertical electron content of the grid at the current moment minus the total vertical electron content of the grid in the previous hour, and then divided by the time interval of 1 hour.

[0018] Preferably, in step 4, the preset area is set to the area with latitude ≥60°N and ≤60°S within the time window corresponding to the suspected abnormal disturbance signal.

[0019] Preferably, in step 4, the Kriging interpolation algorithm is the ordinary Kriging interpolation method, the interpolation lag distance is set to 110km, the search radius is set to 1100km, and a spherical model is used to fit the experimental semivariogram.

[0020] Preferably, in step 4, if the overall disturbance index of the region is greater than three times the average disturbance index of the vertical total electron content spatial region under normal conditions and the proportion of grids exceeding the threshold exceeds 30%, it is determined to be solar or geomagnetic activity interference; otherwise, it is determined to be a local volcanic eruption.

[0021] Preferably, in step 4, the overall disturbance index of the region is the sum of the absolute values ​​of the rate of change of the total vertical electron content of each grid in the vertical total electron content spatial region, divided by the total number of grids in the region.

[0022] Preferably, in step 4, the percentage of grids exceeding the threshold is the proportion of the number of grids in the vertical total electron content spatial region where the rate of change of vertical total electron content exceeds the dynamic threshold of the abnormal state, to the total number of grids in that region.

[0023] Preferably, in step 1, the volcano-related location and eruption data are obtained from the Smithsonian Institution's Global Volcanic Activity Program database, and the occultation data are obtained from the Fengyun Satellite Remote Sensing Data Service Network and the official website of the COSMIC Data Storage Center.

[0024] Preferably, in step 5, the spherical harmonic function used in steps 1 and 3 is of order 15, and the spherical harmonic coefficients are solved by the least squares method to fit the global grid point VTEC.

[0025] Preferably, in step 2, the dynamic threshold for abnormal states and the dynamic threshold for normal states are the sum of the change rate of each grid and three times its standard deviation under the corresponding states.

[0026] The above algorithm effectively solves the problems of complex data processing, limited coverage, and susceptibility to space weather interference in existing methods, enabling rapid and accurate monitoring of volcanic eruptions.

[0027] This addresses the problems of existing ground-based GNSS-based monitoring methods, such as complex and time-consuming data processing, difficulty in achieving rapid global response, and high false alarm rates due to interference from solar and geomagnetic activity.

[0028] This invention has the following characteristics and beneficial effects:

[0029] The method described in this invention enables real-time monitoring of volcanic eruptions using existing COSMIC and Fengyun satellite ionospheric occultation data. Compared to existing ground-based GNSS-based volcanic eruption monitoring methods, it solves the problems of complex data processing and slow response speed leading to difficulties in large-scale detection, while also achieving global dynamic monitoring of volcanic eruptions. Compared to monitoring methods based on a single satellite data source, it solves the problems of limited coverage and difficulty in achieving both spatiotemporal resolution, while avoiding the drawbacks of high false alarm rates caused by interference from solar and geomagnetic activity. By fusing COSMIC and Fengyun satellite occultation data and using spherical harmonic functions to achieve global 1°×1° gridded interpolation, it fully leverages the global coverage advantage of satellites, significantly improving monitoring coverage and spatial resolution. By constructing background models of normal and abnormal ionospheric states, it provides an accurate benchmark for judging ionospheric disturbances caused by volcanic eruptions; by calculating the VTEC rate of change and using dynamic thresholds, it enables rapid identification of suspected abnormal signals, greatly improving monitoring timeliness. The VTEC model, representing the vertical total electron content, was modeled using the Kriging interpolation algorithm. Multi-indicator verification accurately distinguished between volcanic eruptions and interference from solar and geomagnetic activity, effectively reducing the false alarm rate. Finally, a quarterly model update mechanism continuously adapted to the long-term changes in the ionosphere, further ensuring the stability and reliability of monitoring accuracy and providing efficient and accurate technical support for disaster prevention and mitigation, and satellite navigation and communication. Attached Figure Description

[0030] Figure 1 This is a flowchart of a volcanic eruption monitoring method that integrates COSMIC and Fengyun satellite ionospheric occultation data in an embodiment of the present invention.

[0031] Figure 2 This is a flowchart illustrating the background model for establishing the normal and abnormal states of the global ionospheric VTEC based on spherical harmonic functions in an embodiment of the present invention.

[0032] Figure 3 This is a flowchart illustrating the extraction of suspected volcanic eruption signals in an embodiment of the present invention.

[0033] Figure 4 This is a flowchart illustrating the secondary verification and interference differentiation process for the spatial region of vertical total electron content in an embodiment of the present invention. Detailed Implementation

[0034] The present invention will now be described in detail with reference to specific embodiments. These embodiments will help those skilled in the art to further understand the present invention, but do not limit the invention in any way. It should be noted that, unless otherwise specified, the embodiments and features described in the present invention can be combined with each other.

[0035] A method for monitoring volcanic eruptions that integrates data from the Fengyun satellite and the COSMIC ionospheric occultation, such as Figure 1 As shown, it includes the following steps:

[0036] Step 1: Construct a background model of normal and anomalous global ionospheric states: Divide the globe into 1°×1° grids and mark grids containing active and dormant volcanoes. Collect historical ionospheric occultation data from COSMIC-1 and Fengyun-3C satellites during periods of volcanic eruptions without a VEI (Volcanic Eruption Index) of level 4 or above. After data preprocessing, remove anomalous data. Interpolate the electron density profile using height splines. Numerically integrate the electron density in the 110-750km altitude range to obtain the oblique total electron content. Convert the oblique total electron content to the vertical total electron content using a mapping function. Group the vertical total electron content data by hour.

[0037] Specifically, in this embodiment, by querying the Smithsonian Institution's Global Volcanic Activity Plan, multi-source historical data without periods of volcanic eruptions of VEI level 4 or higher were selected. Dates within 30 days after VEI level 4 or higher eruptions were removed to reduce the impact of volcanic eruptions on the ionosphere. The globe was divided into a 1°×1° grid, and grids containing active and dormant volcanoes were marked, denoted as [grid name missing]. Ionospheric occultation data from the Fengyun-3 FY-3C satellite (2015-2019) were downloaded from the Fengyun Satellite Remote Sensing Data Service Network, and occultation data from the COSMIC-1 satellite (during the same period) were downloaded from the COSMIC Data Storage Center (CDDAC) website. The occultation data mainly includes the time, latitude and longitude, electron density, and corresponding altitude of the occultation events. Valid electron density profiles were selected based on the mean deviation (MD) and noise level factor (Delta) for subsequent accurate VTEC calculation. The mean deviation of valid electron density profiles should be 0-1.5, and the noise level factor should be less than 0.01. The formula for calculating the mean deviation is as follows:

[0038]

[0039] in This represents the total number of samples in the electron density profile. For the first Measured electron density at a certain height The background electron density is obtained by performing a 9-point moving average on the measured profile; the noise level factor is calculated using the following formula:

[0040]

[0041] in This represents the total number of electron density sampling points above 300 km. For distances over 300 km Measured electron density at a certain height The corresponding background electron density, This represents the peak electron density of the F2 layer of the ionosphere.

[0042] For a valid electron density profile, spline interpolation is performed at 1 km height intervals to obtain a continuous electron density profile. The TEC value is calculated using the trapezoidal integration method within the height range of 110-750 km (if the highest height of the profile is below 750 km, the actual highest height is used). Since the position of the occultation point gradually shifts with decreasing altitude, the occultation observation provides an electron density profile along a diagonal line in space. Therefore, the TEC calculated by numerical integration is the diagonal TEC, i.e., STEC. The STEC integrated for the k-th electron density profile is denoted as STEC. k Following the processing method of GPS TEC, this method uses the location coordinates at an altitude of 110 km as the observation point coordinates, and introduces a mapping function to process the STEC within the altitude range of 110-750 km. k Converted to vertical total electron content scatter plot ,in Represent the latitude, longitude, and time of the occultation location, respectively, and calculate the zenith angle at the point of penetration. k ,make ,but It can be represented as:

[0043]

[0044] The calculated Data is saved in groups of "year-month-day-hour".

[0045] Furthermore, for each set of data, a spherical harmonic function of order 15 is used for fitting, and the spherical harmonic coefficients are solved using the least squares algorithm, as shown in the following formula:

[0046]

[0047] In the formula The geographical or geomagnetic latitude of the puncture point. The fixed longitude of the puncture point. Let be the highest order of the spherical function (SH) expansion. For the classical, incompletely normalized Legendre function, and These are unknown spherical harmonic coefficients, i.e., the global or regional VTEC parameters to be determined.

[0048] Finally, by substituting the hourly calculated spherical harmonic coefficients and the daily fixed longitude and geographic latitude of the target grid back into equation (4), the vertical total electron content of the target grid can be fitted, denoted as . in for Time, Grid of value.

[0049] Step 2: Define dates with a solar radiation flux index of not less than 150 sfu and an average geomagnetic index of not less than 4 as ionospheric anomalies, and the rest as normal states. Calculate the standard deviation of the hourly rate of change of the vertical total electron content at each grid point under both states. Set the background value of the corresponding spatiotemporal rate of change of the vertical total electron content and the sum of three times the standard deviation as the dynamic thresholds for normal and anomalous states.

[0050] Specifically, such as Figure 2 As shown, for the daily data filtered from 2015 to 2019, the average F10.7 index and the hourly average Kp index were calculated. If F10.7 ≥ 150sfu and the average Kp index ≥ 4, the day was marked as "abnormal state"; otherwise, it was marked as "normal state". The VTEC change rate for each 1°×1° grid and the hourly rate were calculated for both states. The formula for calculating the VTEC change rate is as follows:

[0051]

[0052] In the formula for Time, Grid VTEC value, The time interval is expressed in hours; its standard deviation is calculated, and the standard deviation for each latitude and longitude point is denoted as σ. 1（,i,j,t) and σ 2(,i,j,t) Set the background value of the VTEC rate of change for the corresponding time and space to +3σ. 1（,i,j,t) and the background value of the VTEC rate of change in the corresponding spacetime +3σ 2(,i,j,t) These represent the dynamic threshold for the normal state and the dynamic threshold for VTEC variation caused by solar and geomagnetic activity, respectively, with thresholds T and T1, respectively. 1(i,j,t) and T 2(i,j,t) The thresholds are grouped and saved for easy retrieval later.

[0053] Step 3: Real-time acquisition of ionospheric occultation data from Fengyun-3D to Fengyun-3G satellites and COSMIC-2 satellite. After preprocessing and inversion, the real-time value of the total vertical electron content of the global grid is obtained. The data update frequency is 60 minutes / time. The rate of change of the total vertical electron content of each grid per hour is calculated. If the rate of change exceeds the dynamic threshold and the grid is a marked volcanic area grid, it is judged as a suspected abnormal disturbance signal and the relevant spatiotemporal information is recorded.

[0054] Specifically, such as Figure 3 As shown, in this embodiment, occultation data from the Fengyun-3 series satellites and the COSMIC-2 satellite are collected in real time via a satellite data receiving terminal. The data update cycle is 60 minutes. Using the same method as in step 1, quality control is performed on the real-time occultation data, removing abnormal contour lines, electron density profile interpolation, TEC integral calculation, and VTEC mapping transformation to obtain the real-time VTEC values ​​for a global 1°×1° grid, denoted as [missing information]. For each grid The rate of change of VTEC within the window is calculated using formula (5), and the mesh is extracted simultaneously. At the current hour's normal dynamic threshold, if and exist The grid is then marked as "suspected volcanic eruption," and its latitude and longitude are recorded. Current time and The value of .

[0055] Step 4: Extract the vertical total electron content data of the preset area within the time window corresponding to the suspected abnormal disturbance signal, and use the Kriging interpolation algorithm to construct the vertical total electron content spatial region of the grid within the preset area. Calculate the grid proportion of the overall disturbance index and the abnormal state dynamic threshold within the time window based on the vertical total electron content spatial region. Determine whether it is solar or geomagnetic activity interference based on the proportion; otherwise, determine whether it is a local volcanic eruption.

[0056] In this embodiment, as Figure 4 As shown, the preset area is set as the region with latitude ≥60°N and ≤60°S within the time window corresponding to the suspected abnormal disturbance signal.

[0057] Furthermore, real-time VTEC data for latitudes ≥60°N and ≤60°S are extracted at the current time. Ordinary Kriging (OrK) is used to construct a 1°×1° grid spatial region for the total electron content in the vertical direction. The VTEC value of each grid within the vertical total electron content spatial region is then obtained. The specific interpolation process is as follows:

[0058] Ordinary Kriging interpolation is based on the assumption of intrinsic stationarity of the data, and the relationship between statistical variables satisfies the following expression:

[0059]

[0060]

[0061]

[0062] In the formula, , , Let these represent the expectation, variance, and covariance operators, respectively; Two-dimensional spatial coordinates; and They are spatial points and The vertical delay of the ionosphere at that location corresponds to the VTEC value; This is the separation distance between spatial points, also known as the lag distance; It is a semi-variogram function used to characterize the spatial autocorrelation structure of data.

[0063] First, based on real-time VTEC observation data of the vertical total electron content spatial region, the experimental semivariogram function under different lag times is calculated, as shown in the following formula:

[0064]

[0065] In the formula, Lag distance The corresponding number of observation point pairs; and They are spatial points and The VTEC observation value at the location. In this invention, the lag distance h is set to 110 km and the search radius is set to 1100 km to balance spatial correlation and computational efficiency.

[0066] A spherical model was used to fit the experimental semivariogram. This model can effectively characterize the spatial correlation features of VTEC in the ionosphere, and its mathematical form is as follows:

[0067]

[0068] In the formula, This is known as the nugget effect, reflecting the discontinuity of the origin caused by microscale changes. The range of influence refers to the extent beyond which there is no significant spatial correlation between observations. This is the sill value, which is the value of the semivariogram when it reaches stability.

[0069] Based on the fitted theoretical semi-variogram, the weighting coefficients are solved using the Lagrange multiplier method. The VTEC estimate for unsampled grid points is obtained by weighted averaging of the observed values, as shown in the following formula:

[0070]

[0071] In the formula, For grid points to be estimated VTEC value at the location; For the first The weighting coefficients of each observation point satisfy the following conditions: And minimize the variance of the estimation error. ; The number of observation points participating in the interpolation; For the first Measured VTEC values ​​at each observation point.

[0072] After reconstructing the VTEC values ​​of a 1°×1° grid within the vertical total electron content spatial region using the above-described ordinary Kriging interpolation process, a secondary verification is performed on the grid in step 3 that is suspected of being a volcanic eruption: the VTEC grid data after high-latitude interpolation is then verified. The rate of change of VTEC for each grid at time t is calculated according to formula (5). Extract the dynamic threshold T for ionospheric anomalies constructed in step 1. 2(i,j,t) Statistical analysis of total vertical electron content within spatial regions The number of grid cells exceeding the threshold is calculated, and the proportion R of the grid cells exceeding the threshold to the total grid cells is calculated. At the same time, the overall perturbation index DI of the region and the average perturbation index DI0 of the vertical total electron content spatial region under normal conditions are calculated, which reflects the average oscillation level of VTEC in the vertical total electron content spatial region. The formula for DI is as follows:

[0073]

[0074] In the formula, A grid set representing the spatial region of total electron content in the vertical direction. This represents the total number of grid cells in the region.

[0075] If DI > 3DI0 and R > 30%, it is determined to be interference from solar or geomagnetic activity, and the suspected volcanic eruption marker is removed; if DI ≤ 3DI0 and R ≤ 30%, it is determined to be ionospheric disturbance caused by local volcanic eruption, and the monitoring result is retained and confirmed.

[0076] Step 5: Add the collected and verified volcanic eruption event data and corresponding ionospheric observation data to the historical database and update the dynamic threshold.

[0077] In this embodiment, the final confirmed volcanic eruption monitoring results, including eruption time, latitude and longitude region, VTEC disturbance intensity ΔVTEC, high latitude verification index DI and R, are presented in a standardized text file and a visualization chart. The chart shows the determined location of the volcanic eruption and the distribution map of the over-threshold grid in the high latitude region.

[0078] Quarterly, field-verified volcanic eruption data are collected from the Smithsonian Global Volcanic Activity Programme database and supplemented to the historical database by combining corresponding COSMIC and Fengyun satellite ionospheric observation data; the above steps are then repeated to update the dynamic threshold T.1(i,j,t) and T 2(i,j,t) This ensures that the model can dynamically adapt to the long-term changes in the ionosphere and continuously optimize monitoring accuracy.

[0079] The foregoing has shown and described the basic principles, main features, and advantages of the present invention. Those skilled in the art should understand that the present invention is not limited to the above embodiments. The embodiments and descriptions in the specification are merely preferred examples and are not intended to limit the invention. Various changes and modifications can be made to the invention without departing from its spirit and scope, and all such changes and modifications fall within the scope of the present invention as claimed. The scope of protection of the present invention is defined by the appended claims and their equivalents.

Claims

1. A method for monitoring volcanic eruptions by integrating data from the Fengyun satellite and the COSMIC ionospheric occultation, characterized in that, Includes the following steps: Step 1: Grid the globe and mark the grids containing active and dormant volcanoes. Collect historical ionospheric occultation data from COSMIC-1 and Fengyun-3C satellites during periods without VEI4 or higher volcanic eruptions. After preprocessing the data, remove outliers. Fit the preprocessed data with a spherical harmonic function to invert and obtain the total vertical electron content of each grid point. Step 2: Define abnormal and normal states based on the solar radiation flux index under the grid, and obtain the dynamic thresholds for abnormal and normal states by calculating the standard deviation of the hourly change rate of the vertical total electron content at each grid point under the two states. Step 3: Real-time acquisition of ionospheric occultation data from Fengyun-3D to Fengyun-3G satellites and COSMIC-2 satellite. After preprocessing and spherical harmonic function fitting, the real-time value of the vertical total electron content of the global grid is obtained by inversion. The data update frequency is 60 minutes / time. The rate of change of the vertical total electron content of each grid per hour is calculated. If the rate of change exceeds the dynamic threshold of the normal state and the grid is a marked volcanic area grid, it is judged as a suspected abnormal disturbance signal and the relevant spatiotemporal information is recorded. Step 4: Extract the vertical total electron content data of a preset area within the time window corresponding to the suspected abnormal disturbance signal. Use the Kriging interpolation algorithm to construct the vertical total electron content spatial region of the grid within the preset area. Calculate the overall disturbance index and the proportion of grids exceeding the threshold within the time window based on the vertical total electron content spatial region. Determine whether the disturbance is caused by solar or geomagnetic activity based on the overall disturbance index and the proportion of grids exceeding the threshold; otherwise, determine whether it is a local volcanic eruption. The proportion of grids exceeding the threshold is the percentage of the total number of grids in the vertical total electron content spatial region whose rate of change in vertical total electron content exceeds the dynamic threshold of the abnormal state, relative to the total number of grids in the region. Step 5: Add the collected and verified volcanic eruption event data and corresponding ionospheric observation data to the historical database and update the dynamic threshold.

2. The method according to claim 1, characterized in that, In step 1, during the gridding process, the globe is divided into 1°×1° grids and grids containing active and dormant volcanoes are marked.

3. The method according to claim 2, characterized in that, In step 1, the preprocessing includes determining the quality of the electron density profile based on the average deviation and noise level factor and removing abnormal data, wherein the average deviation of the effective electron density profile is 0-1.5 and the noise level factor is less than 0.

01.

4. The method according to claim 2, characterized in that, In step 1, the method for calculating the vertical total electron content is as follows: the electron density profile in the preprocessed ionospheric occultation data is interpolated using height splines, and the electron density in the 110-750km altitude range is numerically integrated to obtain the oblique total electron content. The oblique total electron content is then converted into vertical total electron content scatter points using a mapping function. The vertical total electron content scatter points are then fitted with a spherical harmonic function to obtain a 1°×1° grid of vertical total electron content.

5. The method according to claim 4, characterized in that, The mapping function is the cosine of the oblique total electron content multiplied by the zenith angle at the puncture point.

6. The method according to claim 1, characterized in that, In step 2, dates with a solar radiation flux index of not less than 150 sfu and an average geomagnetic index of not less than 4 are defined as ionospheric anomalous states, while the rest are normal states.

7. The method according to claim 1, characterized in that, In step 2, the rate of change of the total vertical electron content per hour is the total vertical electron content of the grid at the current moment minus the total vertical electron content of the grid in the previous hour, and then divided by the time interval of 1 hour.

8. The method according to claim 1, characterized in that, In step 4, the preset area is set as the area with latitude ≥60°N and ≤60°S within the time window corresponding to the suspected abnormal disturbance signal.

9. The method according to claim 1, characterized in that, In step 4, the Kriging interpolation algorithm is the ordinary Kriging interpolation method, the interpolation lag distance is set to 110km, the search radius is set to 1100km, and a spherical model is used to fit the experimental semivariogram.

10. The method according to claim 1, characterized in that, In step 4, if the overall disturbance index of the region is greater than three times the average disturbance index of the vertical total electron content spatial region under normal conditions and the proportion of grids exceeding the threshold exceeds 30%, it is determined to be solar or geomagnetic activity interference; otherwise, it is determined to be a local volcanic eruption.

11. The method according to claim 1, characterized in that, In step 4, the overall regional disturbance index is the sum of the absolute values ​​of the rate of change of the total vertical electron content of each grid in the vertical total electron content spatial region, divided by the total number of grids in that region.

12. The method according to claim 1, characterized in that, In step 1, the volcano-related location and eruption data are obtained from the Smithsonian Institution's Global Volcanic Activity Program database, and the occultation data are obtained from the Fengyun Satellite Remote Sensing Data Service Network and the official website of the COSMIC Data Storage Center.

13. The method according to claim 1, characterized in that, In steps 1 and 3, the spherical harmonic function is of order 15, and the spherical harmonic coefficients are solved by the least squares method, thereby fitting the global grid point VTEC.

14. The method according to claim 1, characterized in that, In step 2, the dynamic threshold for abnormal states and the dynamic threshold for normal states are the sum of the change rate of each grid and three times its standard deviation under the corresponding states.

Images