Method and system for determining abnormal disturbance of volcanic ionized layer, storage medium and equipment

Abstract

The embodiment of the invention discloses a method for determining abnormal disturbance of a volcano ionized layer. The method comprises the following steps: acquiring ground-based GNSS (Global Navigation Satellite System) original observation data; preprocessing the ground-based GNSS original observation data to obtain an ionized layer TEC data set; according to the variance of the wavelet energy spectrum measurement time sequence of the ionized layer TEC data set in the time domain and the frequency domain, ionized layer TEC data highly related to the space weather environment is removed, and an optimal ionized layer TEC data set not affected by the space weather environment is obtained; an LSTM algorithm is improved through a Prophet algorithm, and an upper boundary line and a lower boundary line of volcano ionized layer TEC abnormal disturbance are determined based on the optimal ionized layer TEC data set; and determining whether the volcano ionosphere has abnormal disturbance or not according to a comparison relationship between the TEC observation data of the ionosphere and the upper boundary line and the lower boundary line of the abnormal disturbance of the TEC of the volcano ionosphere. According to the method, data processing is carried out, space weather interference is removed, an LSTM algorithm is improved by using a Prophet algorithm to realize abnormal disturbance identification, and volcanic eruption can be efficiently monitored.

Description

Technical Field

[0001] The present invention relates to the technical field of ionospheric abnormal disturbance detection, and in particular to a method, system, storage medium and device for determining volcanic ionospheric abnormal disturbance. Background Art

[0002] Volcanic eruption is a geological process in which magma, volcanic gases and debris from the earth's interior are ejected outward through openings in the earth's crust. This process not only shapes new geological structures, such as spectacular volcanic cones and vast lava plains, but also has a profound impact on the environment and human society. Volcanic activity can trigger climate change, destroy ecosystems, and directly threaten human safety; therefore, monitoring volcanic activity and accurately predicting volcanic eruptions become the key to reducing their potential hazards.

[0003] Volcanic eruptions are not limited to surface activities, their strong energy release can also disturb the Earth's space atmosphere. Specifically, the impact of volcanic eruptions can trigger mesosphere gravity waves (mesospheric airglow waves), which further trigger volcanic ionospheric disturbances (CVIDs), which usually occur within 10 to 45 minutes after the volcanic eruption and present in a quasi-periodic form. These ionospheric disturbances have a significant impact on the total electron content (TEC) in the ionosphere, providing an important basis for the monitoring and prediction of volcanic activities.

[0004] TEC in the ionosphere is a key parameter reflecting the state of the ionosphere, and its changes are closely related to volcanic activity. The ionospheric disturbances generated during and after a volcanic eruption can, to a certain extent, reflect the characteristics and intensity of the volcanic eruption. More importantly, the abnormal disturbances in the ionosphere before a volcanic eruption often indicate the approach of volcanic activity, which is of great significance for the early warning of volcanic eruptions. Summary of the invention

[0005] Based on this, it is necessary to propose a method for determining abnormal disturbances in the volcanic ionosphere to address the above problems.

[0006] A method for determining abnormal disturbance of volcanic ionosphere, the method comprising the following steps:

[0007] Obtain ground-based GNSS raw observation data;

[0008] Preprocessing the ground-based GNSS raw observation data to obtain an ionospheric TEC data set;

[0009] According to the variance of the wavelet energy spectrum of the ionospheric TEC data set in the time domain and the frequency domain, the ionospheric TEC data with high correlation with the space weather environment are removed to obtain the optimal ionospheric TEC data set not affected by the space weather environment;

[0010] The LSTM algorithm is improved by the Prophet algorithm, and the upper and lower limits of the abnormal disturbance of the volcanic ionosphere TEC are determined based on the optimal ionosphere TEC data set;

[0011] According to the comparison relationship between the ionospheric TEC observation data and the upper limit and the lower limit of the abnormal disturbance of the volcanic ionosphere TEC, it is determined whether the volcanic ionosphere has abnormal disturbance.

[0012] In the above scheme, the preprocessing of the ground-based GNSS raw observation data to obtain the ionospheric TEC data set specifically includes:

[0013] Based on the GNSS original pseudorange and carrier phase observation equation, the slant ionospheric delay Iono is determined, where Ion0 = STEC + DCB, STEC is the slant TEC, and VTEC is the vertical TEC;

[0014] Determine the differential code bias DCB based on the non-differential non-combined precise point positioning model and the generalized trigonometric series model;

[0015] The oblique ionospheric delay Iono and the differential code bias DCB are used to determine the oblique TEC and vertical TEC of the satellite puncture point of the GNSS site, and to obtain an ionospheric TEC data set.

[0016] In the above scheme, the determination of the oblique ionospheric delay Iono based on the GNSS original pseudorange and carrier phase observation equation specifically includes:

[0017] The GNSS original pseudorange and carrier phase observation equations are:

[0018]

[0019]

[0020] Wherein, the superscript s and subscript r represent the satellite and receiver respectively, and the subscript f represents the L1 or L2 frequency band of the carrier signal; represents the pseudorange observation value, Represents the carrier phase observation value; represents the geometric distance between the satellite and the receiver; c represents the speed of light; dt r and dt s They represent the clock difference between the receiver and the satellite respectively; Trop represents the tropospheric slant delay, represents the ionospheric slant delay; b r,f and They represent the pseudorange hardware delays at the receiver and satellite ends respectively; f represents the carrier phase wavelength, Indicates the carrier phase integer ambiguity; B r,fand Respectively represent the phase hardware delay of the receiver and satellite ends; ε P and ε L denote the sum of multipath effects, observation noise and other unmodeled errors in pseudorange and carrier observations, respectively;

[0021] The GNSS original pseudorange and carrier phase observation equations are solved according to the least squares method to determine the oblique ionospheric delay Iono.

[0022] In the above scheme, the differential code deviation DCB is determined according to the non-differential non-combined precise point positioning model and the generalized trigonometric series model, specifically including:

[0023] Introduce IGS to publish precise ephemeris and clock products and use ionospheric elimination combination to estimate satellite clock error:

[0024]

[0025] To convert STEC to VTEC, the conversion formula is as follows:

[0026] VTEC=MF·STEC

[0027]

[0028] Where R is the radius of the earth; H is the height of the ionosphere piercing point; α is the coefficient; z is the satellite altitude angle;

[0029] The DCB of the GNSS station is estimated using generalized trigonometric series. The function model of generalized trigonometric series is as follows:

[0030]

[0031] In the formula, and λ represent the geographical latitudes of the ionospheric piercing point and the modeling center point, respectively; h represents the function related to the local time at the ionospheric piercing point; nmax, mmax and kmax represent the maximum orders of the polynomial function and the trigonometric series function, respectively; Enm, Ck and Sk represent the model parameters to be estimated of the generalized trigonometric series function, respectively.

[0032] In the above scheme, the variance of the time series in the time domain and the frequency domain is measured according to the wavelet energy spectrum of the ionospheric TEC data set, the ionospheric TEC data with high correlation with the space weather environment is removed, and the optimal ionospheric TEC data set not affected by the space weather environment is obtained, which specifically includes:

[0033] The mother wavelet φ(t) is simply stretched and translated to obtain a family of sub-wavelets:

[0034]

[0035] Among them, a represents the translation parameter that controls the position of the wavelet, b represents the scale parameter that controls the scaling or expansion of the wavelet width, and t represents time;

[0036] By giving the ionospheric TEC time series or solar activity parameter time series, its continuous wavelet variation formula is obtained:

[0037]

[0038] In the formula, φ * (t) represents the complex conjugate function of φ(t), x represents the ionospheric TEC time series or the solar activity parameter time series;

[0039] Determine the variance of the wavelet energy spectrum measure time series of the ionospheric TEC dataset in the time and frequency domains:

[0040] |(WPS) x (a,b)|=|W x (a,b)| 2

[0041]

[0042] Among them, W x (a, b) and W y (a,b) represent the wavelet transform of x(t) and y(t), respectively. W y The complex conjugate of (a,b).

[0043] In the above scheme, the LSTM algorithm is improved by the Prophet algorithm, and the upper and lower limits of the abnormal disturbance of the volcanic ionosphere TEC are determined based on the optimal ionosphere TEC data set. Before that, the method also includes: obtaining the TEC prediction values ​​of the Prophet algorithm and the LSTM algorithm:

[0044] y(t)=g(t)+s(t)+h(t)+ε t

[0045] In the formula, g(t) represents the trend term; s(t) represents the period term; h(t) represents holidays, which indicates the impact of holidays that may occur irregularly within one or more days; ε t represents the error term; y(t) represents the TEC prediction value of the Prophet algorithm;

[0046] Related parameters for constructing LSTM neural network:

[0047]

[0048]

[0049]

[0050]

[0051]

[0052] In the formula, C t-1 represents the cell state at time t-1, h t-1 represents the output at the previous moment, x t Indicates the current input, h t represents the current output, i.e. the TEC prediction value of the LSTM algorithm, C t Represents the cell state at the current time t.

[0053] In the above scheme, the LSTM algorithm is improved by the Prophet algorithm, and the upper and lower limits of the abnormal disturbance of the volcanic ionosphere TEC are determined based on the optimal ionosphere TEC data set, specifically including:

[0054] The predicted values ​​of the Prophet algorithm are respectively given and the LSTM model predicts the value The initial weights ω1 and ω2 of ionospheric TEC at time i;

[0055] By analyzing the error accuracy with the measured ionospheric TEC, the objective function is determined:

[0056]

[0057] Where N represents the total number of detection periods; yi represents the measured ionospheric TEC; represents the TEC forecast value of the Prophet algorithm at time i, Represents the TEC prediction value of the LSTM algorithm at time i;

[0058] Preset constraints: ω1+ω2=1, ω1≥0, ω2≥0;

[0059] The optimal weights ω1 and ω2 are determined by the least squares method, and the upper and lower limits of the volcanic ionospheric TEC anomaly disturbance are determined.

[0060] In the above scheme, the least square method is used to determine the optimal weights ω1 and ω2, and to determine the upper and lower limits of the abnormal disturbance of the volcanic ionosphere TEC, specifically including:

[0061] Error=TEC actual -TEC predicted

[0062] UB=Mean(Error)+2×Std(Error)

[0063] LB=Mean(Error)-2×Std(Error)

[0064] In the formula, TEC actual Represents the measured ionospheric TEC; TEC predicted represents the ionospheric TEC predicted by the new algorithm; Mean represents the average value; Std represents the standard deviation; UB represents the upper limit of the ionospheric TEC anomaly; LB represents the lower limit of the ionospheric TEC anomaly.

[0065] The present application also proposes a volcanic ionosphere abnormal disturbance detection system, the system comprising: a data acquisition unit, a judgment condition preset unit and a determination unit;

[0066] The data acquisition unit is used to acquire ground-based GNSS raw observation data; pre-process the ground-based GNSS raw observation data to acquire an ionospheric TEC data set; measure the variance of the time series in the time domain and the frequency domain according to the wavelet energy spectrum of the ionospheric TEC data set, remove the ionospheric TEC data with high correlation with the space weather environment, and acquire an optimal ionospheric TEC data set that is not affected by the space weather environment;

[0067] The judgment condition preset unit is used to improve the LSTM algorithm by the Prophet algorithm, and determine the upper limit and the lower limit of the abnormal disturbance of the volcanic ionosphere TEC based on the optimal ionosphere TEC data set;

[0068] The determination unit is used to determine whether the volcanic ionosphere has abnormal disturbance according to the comparison relationship between the ionosphere TEC data set and the upper limit and the lower limit of the abnormal disturbance of the volcanic ionosphere TEC.

[0069] The present application also proposes a readable storage medium storing a computer program. When the computer program is executed by a processor, the processor executes the following steps:

[0070] Obtain ground-based GNSS raw observation data;

[0071] Preprocessing the ground-based GNSS raw observation data to obtain an ionospheric TEC data set;

[0072] According to the variance of the wavelet energy spectrum of the ionospheric TEC data set in the time domain and the frequency domain, the ionospheric TEC data with high correlation with the space weather environment are removed to obtain the optimal ionospheric TEC data set not affected by the space weather environment;

[0073] The LSTM algorithm is improved by the Prophet algorithm, and the upper and lower limits of the abnormal disturbance of the volcanic ionosphere TEC are determined based on the optimal ionosphere TEC data set;

[0074] According to the comparison relationship between the ionospheric TEC observation data and the upper limit and the lower limit of the abnormal disturbance of the volcanic ionosphere TEC, it is determined whether the volcanic ionosphere has abnormal disturbance.

[0075] The present application also proposes a computer device, comprising a memory and a processor, wherein the memory stores a computer program, and the processor executes the following steps of the computer program:

[0076] Obtain ground-based GNSS raw observation data;

[0077] Preprocessing the ground-based GNSS raw observation data to obtain an ionospheric TEC data set;

[0078] According to the variance of the wavelet energy spectrum of the ionospheric TEC data set in the time domain and the frequency domain, the ionospheric TEC data with high correlation with the space weather environment are removed to obtain the optimal ionospheric TEC data set not affected by the space weather environment;

[0079] The LSTM algorithm is improved by the Prophet algorithm, and the upper and lower limits of the abnormal disturbance of the volcanic ionosphere TEC are determined based on the optimal ionosphere TEC data set;

[0080] According to the comparison relationship between the ionospheric TEC observation data and the upper limit and the lower limit of the abnormal disturbance of the volcanic ionosphere TEC, it is determined whether the volcanic ionosphere has abnormal disturbance.

[0081] The embodiment of the present invention has the following beneficial effects: firstly, the original ground-based GNSS observation data is obtained; the original ground-based GNSS observation data is preprocessed to obtain the ionospheric TEC data set; the variance of the time series in the time domain and the frequency domain of the wavelet energy spectrum measurement of the ionospheric TEC data set is removed, and the ionospheric TEC data with high correlation with the space weather environment is obtained, and the optimal ionospheric TEC data set not affected by the space weather environment is obtained; the LSTM algorithm is improved by the Prophet algorithm, and the upper and lower limits of the abnormal disturbance of the volcanic ionosphere TEC are determined based on the optimal ionosphere TEC data set; according to the comparison relationship between the ionospheric TEC observation data and the upper and lower limits of the abnormal disturbance of the volcanic ionosphere TEC, it is determined whether the volcanic ionosphere has abnormal disturbance. The present invention can efficiently monitor volcanic eruptions by processing data and removing space weather interference, using the Prophet algorithm combined with the LSTM algorithm to realize abnormal disturbance identification. BRIEF DESCRIPTION OF THE DRAWINGS

[0082] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative work.

[0083] in:

[0084] Figure 1 It is a schematic flow chart of a method for determining abnormal disturbance of volcanic ionosphere in one embodiment;

[0085] Figure 2 A schematic diagram of the wavelet coherence spectrum of the TEC time series and the space weather parameters of the TONG station in an embodiment;

[0086] Figure 3 A schematic diagram of the wavelet coherence spectrum of the TEC time series and the space weather parameters of the TONG station in an embodiment;

[0087] Figure 4 A schematic diagram of a TEC error sequence change of a SAMO station in one embodiment;

[0088] Figure 5 A schematic diagram of TEC sequence variation of a SAMO station in one embodiment. DETAILED DESCRIPTION

[0089] The technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, rather than all the embodiments; based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative work are within the scope of protection of the present invention.

[0090] In the following description, a large number of specific details are given in order to provide a more thorough understanding of the present invention; however, it is obvious to those skilled in the art that the present invention can be implemented without one or more of these details; in other examples, in order to avoid confusion with the present invention, some technical features known in the art are not described, and it should be understood that the present invention can be implemented in different forms and should not be interpreted as limited to the embodiments presented here; on the contrary, providing these embodiments will make the disclosure thorough and complete and fully convey the scope of the present invention to those skilled in the art.

[0091] The purpose of the terms used herein is only to describe specific embodiments and is not intended to be limiting of the present invention. When used herein, the singular forms "a", "an" and "said / the" are also intended to include the plural forms, unless the context clearly indicates otherwise, and it should also be understood that the terms "consisting of" and / or "comprising", when used in this specification, determine the presence of the features, integers, steps, operations, elements and / or parts, but do not exclude the presence or addition of one or more other features, integers, steps, operations, elements, parts and / or groups. When used herein, the term "and / or" includes any and all combinations of the relevant listed items.

[0092] To facilitate understanding, the relevant terms involved in this application are first introduced below.

[0093] (1) Ionospheric TEC (Total Electron Content) is a key parameter in ionospheric physics research. It describes the total number of all electrons in a column per unit area in the ionosphere. It includes two main forms: VTEC (Vertical Total Electron Content) and STEC (Slant Total Electron Content).

[0094] (2) The Prophet algorithm is a tool for time series forecasting that can efficiently process and forecast time series data with complex trends and seasonal changes. The algorithm achieves rapid iterative optimization of time series models by combining multiple model components with strong interpretability.

[0095] In the existing technology, with the rapid development of space detection technologies such as the Global Navigation Satellite System (GNSS) and occultation, the detection and research of ionospheric anomalies has become a hot topic in the scientific community at home and abroad. However, although these technologies provide rich observation data, there is still a lack of a systematic and efficient method in the field of detection of abnormal ionospheric disturbances before volcanic eruptions. Most of the current methods focus on the analysis of ionospheric disturbances after volcanic eruptions, while the early warning research before volcanic eruptions is relatively weak.

[0096] Therefore, developing a method that can accurately detect abnormal ionospheric disturbances before volcanic eruptions is of great significance for improving the early warning capability of volcanic eruptions and reducing their potential hazards. The method proposed in this paper is based on this background and aims to achieve effective detection and early warning of abnormal ionospheric disturbances before volcanic eruptions by comprehensively utilizing ground-based GNSS data, space weather data and advanced data processing and analysis techniques.

[0097] In order to thoroughly understand the present invention, a detailed structure will be presented in the following description to illustrate the technical solution proposed by the present invention; optional embodiments of the present invention are described in detail as follows, but in addition to these detailed descriptions, the present invention may also have other implementation methods.

[0098] like Figure 1 As shown, in one embodiment, a method for determining abnormal disturbance of volcanic ionosphere is provided, and the method for determining abnormal disturbance of volcanic ionosphere includes steps S101 to S105, which are described in detail as follows:

[0099] S101, obtaining ground-based GNSS original observation data;

[0100] Due to its global coverage and high precision, the use of ground-based GNSS (Global Navigation Satellite System) as a data source provides a reliable data basis for the monitoring of ionospheric TEC (total electron content).

[0101] S102, preprocessing the ground-based GNSS original observation data to obtain an ionospheric TEC data set;

[0102] By preprocessing the original observation data, noise and outliers can be effectively removed, data quality can be improved, and a more accurate ionospheric TEC dataset can be obtained.

[0103] In some embodiments, preprocessing the ground-based GNSS raw observation data to obtain an ionospheric TEC dataset specifically includes:

[0104] Based on the GNSS original pseudorange and carrier phase observation equation, the slant ionospheric delay Iono is determined, where Ion0 = STEC + DCB, STEC is the slant TEC, and VTEC is the vertical TEC;

[0105] Determine the differential code bias DCB based on the non-differential non-combined precise point positioning model and the generalized trigonometric series model;

[0106] The oblique TEC and vertical TEC of the satellite puncture point of the GNSS station are determined according to the oblique ionospheric delay Iono and the differential code bias DCB, and the ionospheric TEC dataset is obtained.

[0107] Preferably, the GNSS original pseudorange and carrier phase observation data (RINEX format file) are used as the data source, the precise ephemeris and clock error products use the post-precision ephemeris and clock error released by the international GNSS service, the differential code bias (DCB) uses the DCB product released by the Chinese Academy of Sciences (CAS), and the TEQC data processing software is used to perform data format standardization and quality inspection analysis on the original measurement data, and then the non-difference non-combined precise single-point positioning method and the generalized trigonometric series method are used to solve and obtain the slant TEC (Slant TEC, STEC) and vertical TEC (Vertical TEC, VTEC) data of the satellite puncture point of the GNSS site.

[0108] In some embodiments, based on the GNSS original pseudorange and the carrier phase observation equation, determining the oblique ionospheric delay Iono specifically includes:

[0109] The GNSS original pseudorange and carrier phase observation equations are:

[0110]

[0111]

[0112] Wherein, the superscript s and subscript r represent the satellite and receiver respectively, and the subscript f represents the L1 or L2 frequency band of the carrier signal; represents the pseudorange observation value, Represents the carrier phase observation value; represents the geometric distance between the satellite and the receiver; c represents the speed of light; dt r and dt s They represent the clock difference between the receiver and the satellite respectively; Trop represents the tropospheric slant delay, represents the ionospheric slant delay; b r,f and They represent the pseudorange hardware delays at the receiver and satellite ends respectively; f represents the carrier phase wavelength, Indicates the carrier phase integer ambiguity; B r,f and Respectively represent the phase hardware delay of the receiver and satellite ends; ε P and ε L denote the sum of multipath effects, observation noise and other unmodeled errors in pseudorange and carrier observations, respectively;

[0113] The GNSS original pseudorange and carrier phase observation equations are solved according to the least squares method to determine the slant ionospheric delay Iono.

[0114] Preferably, the above formulas (1) and (2) are processed to determine the corresponding linearized equations:

[0115]

[0116]

[0117] Where u represents the coefficient matrix of the unknown number after linearization, and x represents the receiver position parameter.

[0118] Specifically, the ionospheric slant delay in the above GNSS raw pseudorange and carrier phase observation equations is That is, the oblique ionospheric delay Iono.

[0119] In some embodiments, determining the differential code bias DCB according to the non-differential non-combined precise point positioning model and the generalized trigonometric series model specifically includes:

[0120] Introduce IGS to publish precise ephemeris and clock products and use ionospheric elimination combination to estimate satellite clock error:

[0121]

[0122] Among them, IGS uses the ionospheric elimination combination to estimate the satellite clock error, and the final clock error product absorbs the satellite pseudorange hardware delay.

[0123] Substituting equation (5) into equation (3) and equation (4), the final non-difference non-combined precise point positioning model can be obtained. The slant ionospheric delay (Iono) estimated by this method includes the receiver DCB, that is, Ion0 = STEC + DCB, so STEC is converted to VTEC, and the conversion formula is as follows:

[0124] VTEC=MF·STEC (6)

[0125]

[0126] Where R is the radius of the earth, which is 6371 km; H is the height of the ionospheric piercing point, which is 450 km, α = 0.9782; z is the satellite altitude angle;

[0127] The DCB of the GNSS station is estimated using generalized trigonometric series. The function model of generalized trigonometric series is as follows:

[0128]

[0129] In the formula, and λ represent the geographical latitudes of the ionospheric piercing point and the modeling center point, respectively; h represents the function related to the local time at the ionospheric piercing point; nmax, mmax and kmax represent the maximum orders of the polynomial function and the trigonometric series function, respectively; Enm, Ck and Sk represent the model parameters to be estimated of the generalized trigonometric series function, respectively.

[0130] Specifically, h = 2π(t-14) / 24, where t is the local time.

[0131] Furthermore, based on the parameters related to the space weather environment before the volcanic eruption, including the 10.7 cm radiation flux (F10.7), extreme ultraviolet radiation (EUV), solar wind speed (Vsw) in the solar activity parameters, the Bz component in the interplanetary magnetic field IMF, the Kp index, Dst index, Ap index, AE index in the geomagnetic activity index, etc., combined with the solved ionospheric TEC, the continuous wavelet transform (CWT) is used to accurately analyze the space weather conditions before the volcanic eruption and eliminate the interference caused by the space weather environment.

[0132] S103, measuring the variance of the time series in the time domain and the frequency domain according to the wavelet energy spectrum of the ionospheric TEC data set, removing the ionospheric TEC data with high correlation with the space weather environment, and obtaining the optimal ionospheric TEC data set that is not affected by the space weather environment;

[0133] By utilizing the characteristics of wavelet analysis in the time and frequency domains, combined with variance analysis, we can effectively identify and remove ionospheric TEC data that are highly correlated with space weather environments (such as solar wind, geomagnetic storms, etc.). This step significantly improves the accuracy and reliability of subsequent analysis and avoids the misjudgment of abnormal volcanic ionosphere disturbances caused by external factors.

[0134] In some embodiments, the variance of the time series of the wavelet energy spectrum of the ionospheric TEC data set in the time domain and the frequency domain is measured, and the ionospheric TEC data with high correlation with the space weather environment is removed to obtain the optimal ionospheric TEC data set not affected by the space weather environment, specifically including:

[0135] The mother wavelet φ(t) is simply stretched and translated to obtain a family of sub-wavelets:

[0136]

[0137] Among them, a represents the translation parameter that controls the position of the wavelet, b represents the scale parameter that controls the scaling or expansion of the wavelet width, and t represents time;

[0138] By giving the ionospheric TEC time series or solar activity parameter time series, its continuous wavelet variation formula is obtained:

[0139]

[0140] In the formula, φ * (t) represents the complex conjugate function of φ(t), x represents the ionospheric TEC time series or the solar activity parameter time series;

[0141] Among them, the position of the wavelet in the time domain is given by a, and the position in the frequency domain is given by b. The wavelet transform maps the original sequence into a function of a and b, and obtains the time and frequency information at the same time.

[0142] Determine the variance of the wavelet energy spectrum measure time series of the ionospheric TEC dataset in the time and frequency domains:

[0143] |(WPS) x (a,b)|=|W x (a,b)| 2 (11)

[0144]

[0145] Among them, W x (a, b) and W y (a,b) represent the wavelet transform of x(t) and y(t), respectively. W y The complex conjugate of (a,b).

[0146] Preferably, in the application, to quantify the relationship between TEC and related parameters, the cross wavelet transform of two time series x(t) and y(t) is defined as shown in equation (4), where W x (a, b) and W y (a,b) are the wavelet transforms of x(t) and y(t), respectively. It is W y The complex conjugate of (a,b):

[0147] like Figure 2 , Figure 3 The figure shows an example of using wavelet transform to eliminate interference from the Sun-Earth space environment in detecting ionospheric anomalies before the Tonga volcano eruption in 2022.

[0148] The correlation between TEC and six solar-terrestrial space environment parameters over TONG station from January 1 to 16 was analyzed by applying cross wavelet transform and wavelet coherence spectrum. Figure 2 and Figure 3 As shown; Figure 2The color bar on the right side of the middle represents the cross-wavelet power spectrum density, and the direction of the arrow represents the phase relationship between the two: rightward indicates that the two sequences are in phase, leftward indicates that they are in opposite phases, vertically downward indicates that the previous sequence changes 1 / 4 cycle ahead of the next sequence, and vertically upward indicates that the next sequence changes 1 / 4 cycle ahead of the previous sequence; Figure 3 The wavelet coherence spectra of TEC and other parameters are given. The color bar on the right is the wavelet coherence value, which represents the strength of the coherence.

[0149] like Figure 2 The figure shows the cross wavelet transform of the TEC time series and space weather parameters of the TONG station from January 1 to 16, 2022. The area enclosed by the thick black solid line passes the standard red noise test at a 95% confidence level, indicating the significance of the cycle; the conical area below the thin black solid line is the wavelet cone of influence (COI) area, which is the area with a greater impact of the edge effect of the wavelet transform data. The thick red solid line indicates the time of the volcanic eruption, as shown in the figure. Figure 3 Shown is a schematic diagram of the wavelet coherence spectrum of the TEC time series and space weather parameters at TONG station from January 1 to 16, 2022.

[0150] In summary, by comparing the resonance period, phase relationship and coherence of TEC and other parameters, the space environment interference on the 6th, 8th, 9th, 14th and 15th can be effectively excluded.

[0151] S104. Improve the LSTM algorithm through the Prophet algorithm and determine the upper and lower limits of the abnormal disturbance of the volcanic ionosphere TEC based on the optimal ionosphere TEC data set;

[0152] The Prophet algorithm is used to improve the LSTM (Long Short-Term Memory Network) algorithm, combining the accuracy of Prophet in time series prediction and the advantages of LSTM in dealing with long-term dependencies. This improvement enables the model to more accurately capture the nonlinear changes and long-term trends in ionospheric TEC data, thereby more accurately determining the upper and lower boundaries of volcanic ionospheric TEC anomaly disturbances.

[0153] In some embodiments, the LSTM algorithm is improved by the Prophet algorithm, and the upper and lower limits of the abnormal disturbance of the volcanic ionosphere TEC are determined based on the optimal ionosphere TEC data set. Before that, the method further includes:

[0154] Get the TEC prediction values ​​of the Prophet algorithm and LSTM algorithm:

[0155] y(t)=g(t)+s(t)+h(t)+ε t (14)

[0156] In the formula, g(t) represents the trend term; s(t) represents the period term; h(t) represents holidays, which indicates the impact of holidays that may occur irregularly within one or more days; ε t represents the error term; y(t) represents the TEC prediction value of the Prophet algorithm;

[0157] Related parameters for constructing LSTM neural network:

[0158]

[0159]

[0160]

[0161]

[0162]

[0163] In the formula, C t-1 represents the cell state at time t-1, h t-1 represents the output at the previous moment, x t Indicates the current input, h t represents the current output, i.e. the TEC prediction value of the LSTM algorithm, C t Represents the cell state at the current time t.

[0164] The Prophet algorithm is a tool for time series prediction, which realizes the rapid iterative optimization of time series models. The algorithm automatically detects trend change points. The LSTM method has a recursive network structure, which is composed of several LSTM modules. The LSTM module includes a memory unit, an input gate, an output gate, and a forget gate, so it has the ability of long-term memory. Integrating the Prophet algorithm into the LSTM network makes the newly constructed method have higher sensitivity and better applicability when detecting ionospheric anomalies.

[0165] In the above, g(t) is used to simulate non-periodic changes in time series values; s(t) represents periodic terms, such as weekly changes and seasonal changes; h(t) represents holidays, which indicates the impact of holidays that may occur irregularly in one or more days; ε t represents the error term, which is assumed to be normally distributed.

[0166] When applying the Prophet algorithm, the principal component analysis method is used to determine the applicability of the three models of the trend term g(t), periodic term s(t) and holiday h(t) in the Prophet algorithm in detecting ionospheric TEC anomalies before volcanic eruptions. At the same time, four types of hyperparameters (trend change degree, seasonal flexibility, holiday flexibility and seasonal multiplication) are adjusted. The initial default values ​​of the four types of hyperparameters are set to 0.05, 10, 10 and additive. The adjustment is made to minimize the RMSE, so as to determine the best detection combination model.

[0167] Preferably, the memory unit in the LSTM network tends to retain useful information in the sequence model, and the most critical forget gate will make the information before the selection forgotten, and its input is the state value c at time t-1 t-1 , the input value x at time t t and the output value h at the previous moment t-1 composition.

[0168] The established LSTM model sequentially includes a bidirectional LSTM layer (units=64), a Dropout layer (rate=0.2), an LSTM layer (units=64), a Dropout layer (rate=0.2) and a Dense layer (units=1). The optimizer is Adam, the loss function is mse, batch=64, epochs=12, validation_split=0.3 are selected when training the model, and the data of the previous 27 days are used to roll-forward the TEC of the next day to obtain the best LSTM model for detecting ionospheric TEC anomalies.

[0169] The Prophet algorithm is well suited for detecting significant ionospheric anomalies, but cannot accurately detect weaker ionospheric anomalies. The LSTM method solves this defect, but has deviations in the stability of detecting anomalies. Therefore, the LSTM algorithm is improved by the Prophet algorithm, and the upper and lower limits of the volcanic ionospheric TEC anomaly disturbance are determined by the optimal ionospheric TEC data set.

[0170] In some embodiments, the LSTM algorithm is improved by the Prophet algorithm, and the upper and lower limits of the abnormal disturbance of the volcanic ionosphere TEC are determined based on the optimal ionosphere TEC data set, specifically including:

[0171] Assign the predicted values ​​of the Prophet algorithm respectively And the LSTM model predicted value The initial weights ω1 and ω2 of ionospheric TEC at time i.

[0172] By analyzing the error accuracy with the measured ionospheric TEC, the objective function is determined:

[0173]

[0174] Where N represents the total number of detection periods; yi represents the measured ionospheric TEC; represents the TEC forecast value of the Prophet algorithm at time i, Represents the TEC prediction value of the LSTM algorithm at time i.

[0175] Preset constraints: ω1+ω2=1, ω1≥0, ω2≥0;

[0176] The optimal weights ω1 and ω2 are determined by the least squares method, and the upper and lower limits of the volcanic ionospheric TEC anomaly disturbance are determined.

[0177] Specifically, y(t) in formula (14) is the predicted value of the Prophet algorithm h in formula (15) t That is, the LSTM model prediction value

[0178] In some embodiments, the optimal weights ω1 and ω2 are determined by the least square method, and the upper and lower limits of the abnormal disturbance of the volcanic ionosphere TEC are determined, specifically including:

[0179] Error=TEC actual -TEC predicted (17)

[0180] UB=Mean(Error)+2×Std(Error) (18)

[0181] LB=Mean(Error)-2×Std(Error) (19)

[0182] In the formula, TEC actual Represents the measured ionospheric TEC; TEC predicted represents the ionospheric TEC predicted by the new algorithm; Mean represents the average value; Std represents the standard deviation; UB represents the upper limit of the ionospheric TEC anomaly; LB represents the lower limit of the ionospheric TEC anomaly.

[0183] S105. Determine whether the volcanic ionosphere has abnormal disturbances based on the comparison between the ionospheric TEC observation data and the upper and lower limits of the abnormal disturbances in the volcanic ionosphere TEC.

[0184] By comparing the ionospheric TEC observation data with the upper and lower limits of the determined volcanic ionospheric TEC abnormal disturbances, it is possible to clearly identify whether volcanic activity has caused abnormal disturbances in the ionosphere. This method is not only scientific and efficient, but can also provide important reference for volcano monitoring and early warning.

[0185] In some embodiments, the values ​​of ω1 and ω2 are set to 0.6 and 0.4 respectively. Figure 4 and Figure 5 The ionospheric anomalies in front of the Tonga volcano in 2022 detected by this method are as follows: Figure 4 Shown is the change of ionospheric TEC over the SAMO station from January 1 to 16, 2022 (the red line indicates the true value of TEC, the blue line indicates the value predicted by the algorithm, and the thick red solid line indicates the time of the volcanic eruption). The thick red solid line marks the specific time of the volcanic eruption. By comparing the red and blue lines, we can clearly see the significant fluctuations in ionospheric TEC before and after the volcanic eruption. Before the volcanic eruption, the algorithm-predicted value (blue line) and the true value (red line) maintained a relatively close trend, but at the time of the volcanic eruption and for a period of time thereafter, the predicted value showed a large deviation, reflecting the strong disturbance of the ionosphere by the volcanic eruption.

[0186] like Figure 5 Shown is: ionospheric TEC anomaly above the SAMO station from January 1 to 16, 2022 (the black line and blue line represent the upper and lower boundaries of the TEC anomaly, respectively, and the thick red solid line represents the time of the volcanic eruption); the thick red solid line also marks the time of the volcanic eruption. It can be seen from the figure that before and after the volcanic eruption, the TEC value quickly crossed the anomaly limit, indicating that volcanic activity has caused significant abnormal disturbances to the ionosphere. By comparing the degree of deviation of the TEC value from the anomaly limit, the specific degree and duration of the impact of the volcanic eruption on the ionosphere can be quantitatively evaluated.

[0187] In summary, before a volcanic eruption, the method proposed in this application can better predict the changing trend of ionospheric TEC. By setting the upper and lower limits of TEC anomalies, the abnormal ionospheric disturbances caused by volcanic eruptions can be effectively identified, proving that the scheme has high sensitivity and accuracy in determining abnormal ionospheric disturbances of volcanoes. By real-time monitoring of changes in ionospheric TEC, signs of volcanic activity can be discovered in advance, providing valuable warning time for relevant departments.

[0188] The present application also proposes a volcanic ionosphere abnormal disturbance detection system, the system comprising: a data acquisition unit, a judgment condition preset unit and a determination unit;

[0189] A data acquisition unit is used to acquire ground-based GNSS raw observation data; pre-process the ground-based GNSS raw observation data to acquire an ionospheric TEC data set; measure the variance of the time series in the time domain and frequency domain according to the wavelet energy spectrum of the ionospheric TEC data set, remove the ionospheric TEC data with high correlation with the space weather environment, and acquire the optimal ionospheric TEC data set that is not affected by the space weather environment;

[0190] The judgment condition preset unit is used to improve the LSTM algorithm through the Prophet algorithm and determine the upper and lower limits of the abnormal disturbance of the volcanic ionosphere TEC based on the optimal ionosphere TEC data set;

[0191] The determination unit is used to determine whether the volcanic ionosphere has abnormal disturbance according to the comparison relationship between the ionosphere TEC data set and the upper limit and the lower limit of the abnormal disturbance of the volcanic ionosphere TEC.

[0192] The present application also proposes a readable storage medium storing a computer program. When the computer program is executed by a processor, the processor executes the following steps:

[0193] Obtain ground-based GNSS raw observation data;

[0194] Preprocess the ground-based GNSS raw observation data to obtain the ionospheric TEC dataset;

[0195] According to the wavelet energy spectrum of the ionospheric TEC data set, the variance of the time series in the time domain and frequency domain is measured, and the ionospheric TEC data with high correlation with the space weather environment is removed to obtain the optimal ionospheric TEC data set that is not affected by the space weather environment.

[0196] The LSTM algorithm is improved by the Prophet algorithm, and the upper and lower limits of the abnormal disturbance of the volcanic ionosphere TEC are determined based on the optimal ionosphere TEC data set;

[0197] Based on the comparison between the ionospheric TEC observation data and the upper and lower limits of the abnormal disturbance of the volcanic ionosphere TEC, it is determined whether the volcanic ionosphere has abnormal disturbance.

[0198] The present application also proposes a computer device, including a memory and a processor, wherein the memory stores a computer program, and the computer program is executed by the processor to perform the following steps:

[0199] Obtain ground-based GNSS raw observation data;

[0200] Preprocess the ground-based GNSS raw observation data to obtain the ionospheric TEC dataset;

[0201] According to the wavelet energy spectrum of the ionospheric TEC data set, the variance of the time series in the time domain and frequency domain is measured, and the ionospheric TEC data with high correlation with the space weather environment is removed to obtain the optimal ionospheric TEC data set that is not affected by the space weather environment.

[0202] The LSTM algorithm is improved by the Prophet algorithm, and the upper and lower limits of the abnormal disturbance of the volcanic ionosphere TEC are determined based on the optimal ionosphere TEC data set;

[0203] Based on the comparison between the ionospheric TEC observation data and the upper and lower limits of the abnormal disturbance of the volcanic ionosphere TEC, it is determined whether the volcanic ionosphere has abnormal disturbance.

[0204] Those skilled in the art can understand that all or part of the processes in the above-mentioned embodiment methods can be completed by instructing the relevant hardware through a computer program, and the program can be stored in a non-volatile computer-readable storage medium. When the program is executed, it can include the processes of the embodiments of the above-mentioned methods. Among them, any reference to memory, storage, database or other media used in the embodiments provided in this application can include non-volatile and / or volatile memory. Non-volatile memory may include read-only memory (ROM), programmable ROM (PROM), electrically programmable ROM (EPROM), electrically erasable programmable ROM (EEPROM) or flash memory. Volatile memory may include random access memory (RAM) or external cache memory. As an illustration and not limitation, RAM is available in many forms, such as static RAM (SRAM), dynamic RAM (DRAM), synchronous DRAM (SDRAM), double data rate SDRAM (DDRSDRAM), enhanced SDRAM (ESDRAM), synchronous link (Synchlink) DRAM (SLDRAM), memory bus (Rambus) direct RAM (RDRAM), direct memory bus dynamic RAM (DRDRAM), and memory bus dynamic RAM (RDRAM).

[0205] The technical features of the above embodiments can be combined arbitrarily. To make the description concise, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this specification.

[0206] The above-mentioned embodiments only express several implementation methods of the present application, and the descriptions thereof are relatively specific and detailed, but they cannot be understood as limiting the scope of the patent of the present application. It should be pointed out that, for ordinary technicians in this field, several modifications and improvements can be made without departing from the concept of the present application, and these all belong to the protection scope of the present application. The above disclosure is only a preferred embodiment of the present invention, and of course it cannot be used to limit the scope of the rights of the present invention. Therefore, equivalent changes made according to the claims of the present invention are still within the scope covered by the present invention.

Claims

1. A method for determining abnormal disturbances in the volcanic ionosphere, characterized in that: The method comprises: Obtain ground-based GNSS raw observation data; Preprocessing the ground-based GNSS raw observation data to obtain an ionospheric TEC data set; According to the variance of the wavelet energy spectrum of the ionospheric TEC data set in the time domain and the frequency domain, the ionospheric TEC data with high correlation with the space weather environment are removed to obtain the optimal ionospheric TEC data set not affected by the space weather environment; The LSTM algorithm is improved by the Prophet algorithm, and the upper and lower limits of the abnormal disturbance of the volcanic ionosphere TEC are determined based on the optimal ionosphere TEC data set; According to the comparison relationship between the ionospheric TEC observation data and the upper limit and the lower limit of the abnormal disturbance of the volcanic ionosphere TEC, it is determined whether the volcanic ionosphere has abnormal disturbance.

2. The method for determining abnormal disturbance of volcanic ionosphere according to claim 1, characterized in that: The preprocessing of the ground-based GNSS raw observation data to obtain the ionospheric TEC data set specifically includes: Based on the GNSS original pseudorange and carrier phase observation equation, the slant ionospheric delay Iono is determined, where Ion0 = STEC + DCB, STEC is the slant TEC, and VTEC is the vertical TEC; Determine the differential code bias DCB based on the non-differential non-combined precise point positioning model and the generalized trigonometric series model; The oblique ionospheric delay Iono and the differential code bias DCB are used to determine the oblique TEC and vertical TEC of the satellite puncture point of the GNSS site, and to obtain an ionospheric TEC data set.

3. The method for determining abnormal disturbance of volcanic ionosphere according to claim 2, characterized in that: The determining of the oblique ionospheric delay Iono based on the GNSS original pseudorange and carrier phase observation equation specifically includes: The GNSS original pseudorange and carrier phase observation equations are: Wherein, the superscript s and subscript r represent the satellite and receiver respectively, and the subscript f represents the L1 or L2 frequency band of the carrier signal; represents the pseudorange observation value, Represents the carrier phase observation value; represents the geometric distance between the satellite and the receiver; c represents the speed of light; dt r and dt s They represent the clock difference between the receiver and the satellite respectively; Trop represents the tropospheric slant delay, represents the ionospheric slant delay; b r,f and They represent the pseudorange hardware delays at the receiver and satellite ends respectively; f represents the carrier phase wavelength, Indicates the carrier phase integer ambiguity; B r,f and Respectively represent the phase hardware delay of the receiver and satellite ends; ε P and ε L denote the sum of multipath effects, observation noise and other unmodeled errors in pseudorange and carrier observations, respectively; The GNSS original pseudorange and carrier phase observation equations are solved according to the least squares method to determine the oblique ionospheric delay Iono.

4. The method for determining abnormal disturbance of volcanic ionosphere according to claim 3, characterized in that: The method of determining the differential code bias DCB according to the non-differential non-combined precise point positioning model and the generalized trigonometric series model specifically includes: Introduce IGS to publish precise ephemeris and clock products and use ionospheric elimination combination to estimate satellite clock error: To convert STEC to VTEC, the conversion formula is as follows: VTEC=MF·STEC Where R is the radius of the earth; H is the height of the ionosphere piercing point; α is the coefficient; z is the satellite altitude angle; The DCB of the GNSS station is estimated using generalized trigonometric series. The function model of generalized trigonometric series is as follows: In the formula, and λ represent the geographical latitudes of the ionospheric piercing point and the modeling center point, respectively; h represents the function related to the local time at the ionospheric piercing point; nmax, mmax and kmax represent the maximum orders of the polynomial function and the trigonometric series function, respectively; Enm, Ck and Sk represent the model parameters to be estimated of the generalized trigonometric series function, respectively.

5. The method for determining abnormal disturbance of volcanic ionosphere according to claim 4, characterized in that: The method measures the variance of the time series in the time domain and the frequency domain according to the wavelet energy spectrum of the ionospheric TEC data set, removes the ionospheric TEC data with high correlation with the space weather environment, and obtains the optimal ionospheric TEC data set not affected by the space weather environment, specifically including: The mother wavelet φ(t) is simply stretched and translated to obtain a family of sub-wavelets: Among them, a represents the translation parameter that controls the position of the wavelet, b represents the scale parameter that controls the scaling or expansion of the wavelet width, and t represents time; By giving the ionospheric TEC time series or solar activity parameter time series, its continuous wavelet variation formula is obtained: In the formula, φ * (t) represents the complex conjugate function of φ(t), x represents the ionospheric TEC time series or the solar activity parameter time series; Determine the variance of the wavelet energy spectrum measure time series of the ionospheric TEC dataset in the time and frequency domains: |(WPS) x (a,b)|=|W x (a,b)| 2 Among them, W x (a, b) and W y (a,b) represent the wavelet transform of x(t) and y(t), respectively. W y The complex conjugate of (a,b).

6. The method for determining abnormal disturbance of volcanic ionosphere according to claim 5, characterized in that: Before improving the LSTM algorithm by the Prophet algorithm and determining the upper and lower limits of the abnormal disturbance of the volcanic ionosphere TEC based on the optimal ionosphere TEC data set, the method further includes: obtaining the TEC prediction values ​​of the Prophet algorithm and the LSTM algorithm: y(t)=g(t)+s(t)+h(t)+ε t Where g(t) represents the trend term; s(t) represents the period term; h(t) represents holidays, which indicates the impact of holidays that may occur irregularly within one or more days; ε t represents the error term; y(t) represents the TEC prediction value of the Prophet algorithm; Related parameters for constructing LSTM neural network: In the formula, C t-1 represents the cell state at time t-1, h t-1 represents the output at the previous moment, x t Indicates the current input, h t represents the current output, i.e. the TEC prediction value of the LSTM algorithm, C t Represents the cell state at the current time t.

7. The method for determining abnormal disturbance of volcanic ionosphere according to claim 6, characterized in that: The method of improving the LSTM algorithm by the Prophet algorithm and determining the upper and lower limits of the abnormal disturbance of the volcanic ionosphere TEC based on the optimal ionosphere TEC data set specifically includes: The predicted values ​​of the Prophet algorithm are respectively given and the LSTM model predicts the value The initial weights ω1 and ω2 of ionospheric TEC at time i; By analyzing the error accuracy with the measured ionospheric TEC, the objective function is determined: Where N represents the total number of detection periods; yi represents the measured ionospheric TEC; represents the TEC forecast value of the Prophet algorithm at time i, Represents the TEC prediction value of the LSTM algorithm at time i; Preset constraints: ω1+ω2=1, ω1≥0, ω2≥0; The optimal weights ω1 and ω2 are determined by the least squares method, and the upper and lower limits of the volcanic ionospheric TEC anomaly disturbance are determined.

8. The method for determining abnormal disturbance of volcanic ionosphere according to claim 7, characterized in that: The method of using the least square method to determine the optimal weights ω1 and ω2 and to determine the upper and lower limits of the abnormal disturbance of the volcanic ionosphere TEC specifically includes: Error=TEC actual -TEC predicted UB=Mean(Error)+2×Std(Error) LB=Mean(Error)-2×Std(Error) In the formula, TEC actual Represents the measured ionospheric TEC; TEC predicted represents the ionospheric TEC predicted by the new algorithm; Mean represents the average value; Std represents the standard deviation; UB represents the upper limit of the ionospheric TEC anomaly; LB represents the lower limit of the ionospheric TEC anomaly.

9. A volcanic ionosphere abnormal disturbance detection system, characterized in that: The system comprises: a data acquisition unit, a judgment condition preset unit and a determination unit; The data acquisition unit is used to acquire ground-based GNSS raw observation data; pre-process the ground-based GNSS raw observation data to acquire an ionospheric TEC data set; measure the variance of the time series in the time domain and the frequency domain according to the wavelet energy spectrum of the ionospheric TEC data set, remove the ionospheric TEC data with high correlation with the space weather environment, and acquire an optimal ionospheric TEC data set that is not affected by the space weather environment; The judgment condition preset unit is used to improve the LSTM algorithm by the Prophet algorithm, and determine the upper limit and the lower limit of the abnormal disturbance of the volcanic ionosphere TEC based on the optimal ionosphere TEC data set; The determination unit is used to determine whether the volcanic ionosphere has abnormal disturbance according to the comparison relationship between the ionosphere TEC data set and the upper limit and the lower limit of the abnormal disturbance of the volcanic ionosphere TEC.

10. A readable storage medium storing a computer program, wherein when the computer program is executed by a processor, the processor is caused to perform the steps of the method according to any one of claims 1 to 8.

11. A computer device comprising a memory and a processor, wherein the memory stores a computer program, and when the computer program is executed by the processor, the processor executes the steps of the method according to any one of claims 1 to 8.

Images