[0001] 本发明涉及电离层异常扰动探测技术领域，尤其涉及火山电离层异常扰动确定方法、系统、存储介质及设备。

[0002] 火山喷发是地球内部岩浆、火山气体及碎屑通过地壳中的开口向外喷出的地质过程，这一过程不仅塑造了新的地质结构，如壮观的火山锥和广阔的熔岩平原，还对环境及人类社会造成了深远的影响，火山活动能够引发气候变迁、破坏生态系统，并直接威胁到人类安全；因此，监测火山活动、准确预测火山喷发成为减轻其潜在危害的关键所在。

[0003] 火山喷发不仅局限于地表活动，其强烈的能量释放还能扰动地球的空间大气环境。具体来说，火山喷发的冲击能够触发中间层重力波(中间层气辉波)，这些波动进一步引发同火山电离层扰动CVID，通常在火山喷发后的10至45分钟内发生，并以准周期性的形式呈现，这些电离层扰动对电离层中的总电子含量(Total Electron Content,TEC)具有显著的影响，为火山活动的监测与预测提供了重要依据。

[0004] 电离层中的TEC是反映电离层状态的关键参数，其变化与火山活动密切相关，火山喷发过程中及之后产生的电离层扰动，能够在一定程度上反映火山喷发的特征和强度，更为关键的是，火山喷发前电离层中出现的异常扰动，往往预示着火山活动的临近，对于火山喷发的预警具有重要意义。

[0089] 下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例；基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

[0090] 在下文的描述中，给出了大量具体的细节以便提供对本发明更为彻底的理解；然而，对于本领域技术人员而言显而易见的是，本发明可以无需一个或多个这些细节而得以实施；在其他的例子中，为了避免与本发明发生混淆，对于本领域公知的一些技术特征未进行描述，应当理解的是，本发明能够以不同形式实施，而不应当解释为局限于这里提出的实施例；相反地，提供这些实施例将使公开彻底和完全，并且将本发明的范围完全地传递给本领域技术人员。

[0091] 在此使用的术语的目的仅在于描述具体实施例并且不作为本发明的限制。在此使用时，单数形式的“一”、“一个”和“所述/该”也意图包括复数形式，除非上下文清楚指出另外的方式，还应明白术语“组成”和/或“包括”，当在该说明书中使用时，确定所述特征、整数、步骤、操作、元件和/或部件的存在，但不排除一个或更多其它的特征、整数、步骤、操作、元件、部件和/或组的存在或添加。在此使用时，术语“和/或”包括相关所列项目的任何及所有组合。

[0092] 为了便于理解，下面先对本申请涉及的相关术语进行介绍。

[0093] (1)电离层TEC，(Total Electron Content，总电子含量)是电离层物理研究中的一个关键参数，它描述了电离层中单位面积柱体内所有电子的数量总和，其包含两种主要的形式：VTEC(Vertical Total Electron Content，垂直总电子含量)和STEC(Slant Total Electron Content，倾斜总电子含量)；

[0094] (2)Prophet算法，是一种用于时间序列预测的工具，它能够高效地处理并预测具有复杂趋势和季节性变化的时间序列数据，该算法通过组合多个可解释性强的模型组件，实现了时间序列模型的快速迭代优化。

[0095] 现有技术中，随着全球导航卫星系统(GNSS)及掩星等空间探测技术的迅猛发展，电离层异常的探测与研究已成为国内外科学界的热点，然而，尽管这些技术提供了丰富的观测数据，但在火山喷发前电离层异常扰动的探测领域，仍缺乏一套系统、高效的方法，目前的方法大多侧重于火山喷发后电离层扰动的分析，而对于火山喷发前的预警研究则显得相对薄弱。

[0096] 因此，开发一种能够准确探测火山喷发前电离层异常扰动的方法，对于提升火山喷发的预警能力、减轻其潜在危害具有重要意义，本文提出的方法正是基于这一背景，旨在通过综合利用地基GNSS数据、空间天气数据以及先进的数据处理与分析技术，实现对火山喷发前电离层异常扰动的有效探测与预警。

[0097] 为了彻底理解本发明，将在下列的描述中提出详细的结构，以便阐释本发明提出的技术方案；本发明的可选实施例详细描述如下，然而除了这些详细描述外，本发明还可以具有其他实施方式。

[0098] 如图1所示，在一个实施例中，提供了一种火山电离层异常扰动确定方法，该火山电离层异常扰动确定方法包括步骤S101至步骤S105，详述如下：

[0099] S101、获取地基GNSS原始观测数据；

[0100] GNSS因其全球覆盖和高精度特性，使用地基GNSS(全球导航卫星系统)作为数据源为电离层TEC(总电子含量)的监测提供了可靠的数据基础。

[0101] S102、对地基GNSS原始观测数据进行预处理，获取电离层TEC数据集；

[0102] 通过对原始观测数据进行预处理，能够有效去除噪声和异常值，提高数据质量，从而获取更准确的电离层TEC数据集。

[0103] 在一些实施例中，对地基GNSS原始观测数据进行预处理，获取电离层TEC数据集，具体包括：

[0104] 基于GNSS原始伪距和载波相位观测方程，确定斜电离层延迟Iono，其中，Ion0＝STEC+DCB，STEC为倾斜TEC，VTEC为垂直TEC；

[0105] 根据非差非组合精密单点定位模型和广义三角级数模型确定差分码偏差DCB；

[0106] 根据斜电离层延迟Iono和差分码偏差DCB确定GNSS站点卫星穿刺点的为倾斜TEC和垂直TEC，获取电离层TEC数据集。

[0107] 优选的，以GNSS原始伪距和载波相位观测数据(RINEX格式文件)为数据源，精密星历及钟差产品使用国际GNSS服务发布的事后精密星历和钟差，差分码偏差(Differential Code Bias，DCB)使用中国科学院(Chinese Academy of Sciences，CAS)发布的DCB产品，应用TEQC数据处理软件对原始测量数据进行数据格式标准化处理及质量检核分析，接着利用非差非组合精密单点定位方法和广义三角级数方法解算获得GNSS站点卫星穿刺点的倾斜TEC(Slant TEC,STEC)和垂直TEC(Vertical TEC,VTEC)数据。

[0108] 在一些实施例中，基于GNSS原始伪距和载波相位观测方程，确定斜电离层延迟Iono，具体包括：

[0109] GNSS原始伪距和载波相位观测方程为：

[0110]

[0111]

[0112] 其中，上标s和下标r分别表示卫星和接收机，下标f表示载波信号的L1或L2频段；表示伪距观测值， 表示载波相位观测值； 表示卫星与接收机之间的几何距离；c表s
示光速；dtr和dt分别表示接收机与卫星的钟差；Trop表示对流层斜延迟， 表示电离层斜延迟；br,f和 分别表示接收机端和卫星端的伪距硬件延迟；λf表示载波相位波长，表示载波相位整周模糊度；Br,f和 分别表示接收机端和卫星端的相位硬件延迟；εP和εL分别表示伪距和载波观测中多路径效应、观测噪声和其他未模型化误差之和；

[0113] 根据最小二乘法对GNSS原始伪距和载波相位观测方程进行求解，确定斜电离层延迟Iono。

[0114] 优选的，对上述公式(1)和(2)进行处理，确定对应的线性化方程：

[0115]

[0116]

[0117] 式中，u表示线性化后未知数的系数阵，x表示接收机位置参数。

[0118] 具体的，上述GNSS原始伪距和载波相位观测方程中的电离层斜延迟 即斜电离层延迟Iono。

[0119] 在一些实施例中，根据非差非组合精密单点定位模型和广义三角级数模型确定差分码偏差DCB，具体包括：

[0120] 引入IGS发布精密星历和钟差产品并采用消电离层组合估计卫星钟差：

[0121]

[0122] 其中，IGS采用消电离层组合估计卫星钟差，最终的钟差产品吸收了卫星伪距硬件延迟。

[0123] 将式(5)带入式(3)和式(4)，可得到最终的非差非组合精密单点定位模型。该方法估计的斜电离层延迟(Iono)包含了接收机DCB，即Ion0＝STEC+DCB，因此将STEC转换为VTEC，转换公式如下：

[0124] VTEC＝MF·STEC (6)

[0125]

[0126] 其中，R表示地球半径，取6371km；H表示电离层穿刺点高度，取值为450km，α＝0.9782；z表示卫星高度角；

[0127] 利用广义三角级数估计GNSS测站的DCB，广义三角级数函数模型如下：

[0128]

[0129] 式中，和λ分别表示电离层穿刺点和建模中心点的地理纬度；h表示电离层穿刺点处与地方时相关的函数；nmax，mmax和kmax分别表示多项式函数及三角级数函数的最大阶次；Enm，Ck和Sk分别表示广义三角级数函数的待估模型参数。

[0130] 具体的，h＝2π(t‑14)/24，t为地方时间。

[0131] 进一步的，基于火山喷发前相关反映空间天气环境的参数，包括太阳活动参数中的10.7厘米辐射通量(F10.7)、极紫外辐射(EUV)、太阳风速度(Vsw)，行星际磁场IMF中的Bz分量，地磁活动指数中的Kp指数、Dst指数、Ap指数、AE指数等，结合解算的电离层TEC，利用连续小波变换(Continuous Wavelet Transform,CWT)精准分析火山喷发前空间天气情况，排除空间天气环境带来的干扰。

[0132] S103、根据电离层TEC数据集的小波能量谱测度时间序列在时域和频域内的方差，去除与空间天气环境相关性高的电离层TEC数据，获取不受空间天气环境影响的最优电离层TEC数据集；

[0133] 利用小波分析在时域和频域内的特性，结合方差分析，能够有效地识别并去除与空间天气环境(如太阳风、地磁暴等)高度相关的电离层TEC数据，这一步骤显著提高了后续分析的准确性和可靠性，避免了外部因素对火山电离层异常扰动的误判。

[0134] 在一些实施例中，根据电离层TEC数据集的小波能量谱测度时间序列在时域和频域内的方差，去除与空间天气环境相关性高的电离层TEC数据，获取不受空间天气环境影响的最优电离层TEC数据集，具体包括：

[0135] 将母小波φ(t)通过简单的伸缩和平移得到子小波族：

[0136]

[0137] 其中，a表示控制小波位置的平移参数，b表示控制小波宽度缩放或扩张的尺度参数，t表示时间；

[0138] 通过给定电离层TEC时间序列或太阳活动参数时间序列，得到其连续小波变化公式：

[0139]

[0140] 式中，φ*(t)表示φ(t)的复共轭函数，x表示电离层TEC时间序列或太阳活动参数时间序列；

[0141] 其中，小波在时域上的位置由a给出，在频域上的位置由b给出，小波变换将原始序列映射为a和b的函数，同时得到时间和频率的信息。

[0142] 确定电离层TEC数据集的小波能量谱测度时间序列在时间和频率域内的方差：

[0143] |(WPS)x(a,b)|＝|Wx(a,b)|2  (11)

[0144]

[0145] 其中，Wx(a,b)和Wy(a,b)分别表示x(t)和y(t)的小波变换， 表示Wy(a,b)的复共轭。

[0146] 优选的，在应用中，量化TEC和相关参数的关系，两个时间序列x(t)和y(t)的交叉小波变换定义见式(4)，其中Wx(a,b)和Wy(a,b)分别是x(t)和y(t)的小波变换， 是Wy(a,b)的复共轭：

[0147] 如图2、图3所示为2022年汤加火山喷发前使用小波变换排除日地空间环境对探测电离层异常的干扰的实例。

[0148] 通过应用了交叉小波变换和小波相干谱分析了1月1日‑16日TONG站上空TEC与六种日地空间环境参数的相关性，如图2和图3所示；图2中右侧颜色条表示交叉小波功率谱密度，箭头方向表示两者的相位关系：向右表示两序列同相位，向左表示反相位，垂直向下表示前一序列提前后一序列1/4个周期变化，垂直向上则表示后一序列比前一序列提前1/4个周期变化；图3给出了TEC与其他参数的小波相干谱，右侧颜色条是小波相干值，表征相干性的强弱。

[0149] 如图2所示为2022年1月1日‑16日TONG站TEC时间序列与空间天气参数的交叉小波变换，其中，黑色粗实线封闭区域通过了95％置信水平的标准红噪声检验，表明周期的显著性；黑色细实线下方锥形区域为小波影响锥(Cone of Influence，COI)区域，为小波变换数据边缘效应影响较大的区域，红色粗实线表示火山喷发时刻，如图3所示为2022年1月1日‑16日TONG站TEC时间序列与空间天气参数的小波相干谱示意图。

[0150] 综上，通过比较TEC与其他参数的共振周期、相位关系及相干性，可以有效排除6、8、9、14和15日的空间环境干扰。

[0151] S104、通过Prophet算法改进LSTM算法，基于最优电离层TEC数据集确定火山电离层TEC异常扰动的上界线、下界线；

[0152] 采用Prophet算法改进LSTM(长短期记忆网络)算法，结合了Prophet在时间序列预测中的准确性和LSTM在处理长期依赖关系上的优势，这种改进使得模型能够更准确地捕捉电离层TEC数据中的非线性变化和长期趋势，从而更精确地确定火山电离层TEC异常扰动的上界线和下界线。

[0153] 在一些实施例中，通过Prophet算法改进LSTM算法，基于最优电离层TEC数据集确定火山电离层TEC异常扰动的上界线、下界线，之前，该方法还包括：

[0154] 获取Prophet算法、LSTM算法的TEC预测值：

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

[0156] 式中，g(t)表示趋势项；s(t)表示周期项；h(t)表示节假日，表示在一天或多天内可能不定期发生的假期的影响；εt表示误差项；y(t)表示Prophet算法的TEC预测值；

[0157] 构建LSTM神经网络的相关参数：

[0158]

[0159]

[0160]

[0161]

[0162]

[0163] 式中，Ct‑1表示t‑1时刻的细胞状态，ht‑1表示前一时刻的输出，xt表示当前输入，ht表示当前输出即LSTM算法的TEC预测值，Ct表示当前t时刻的细胞状态。

[0164] Prophet算法是用于时序预测的工具，其实现了时间序列模型的快速迭代优化，对于趋势变化点，该算法会自动检测，LSTM方法具有递归的网络结构，它由若干个LSTM模块组成，LSTM模块包括记忆单元、输入门、输出门和遗忘门，从而具有长期记忆的能力，将Prophet算法融入LSTM网络中，使新构建的方法探测电离层异常时有更高的灵敏度和更优的适用性。

[0165] 上述，g(t)用于模拟时间序列值中的非周期性变化；s(t)表示周期项，例如一周变化、季节变化；h(t)表示节假日，表示在一天或多天内可能不定期发生的假期的影响；εt代表误差项，假设符合正态分布。

[0166] 在应用Prophet算法时，利用主成分分析法确定Prophet算法中趋势项g(t)、周期项s(t)和节假日h(t)三种模型在探测火山喷发前电离层TEC异常的适用性，同时调整四类超参数(趋势变化程度、季节灵活性、假日灵活性和季节乘法)，四类超参数初始默认值设为0.05、10、10、additive，调整使得RMSE最下，从而确定最佳探测组合模型。

[0167] 优选的，LSTM网络中的记忆单元往往会保留序列模型中的有用信息，而最关键的遗忘门会让选择之前的信息遗忘，其输入是由t‑1时刻的状态值ct‑1，t时刻的输入值xt和上一时刻的输出值ht‑1组成。

[0168] 建立的LSTM模型顺序包括一个双向LSTM层(units＝64)、一个Dropout层(rate＝0.2)、一个LSTM层(units＝64)、一个Dropout层(rate＝0.2)和一个Dense层(units＝1)，优化器选择为Adam，损失函数为mse，训练模型时选择batch＝64，epochs＝12，validation_split＝0.3，利用前27天的数据滚动预测下一天的TEC，获得最佳探测电离层TEC异常的LSTM模型。

[0169] Prophet算法对探测显著的电离层异常有较好地适用性，但无法准确探测较微弱的电离层异常；而LSTM方法解决了这一缺陷，但对探测异常的稳定性有偏差，因此，通过Prophet算法改进LSTM算法，通过最优电离层TEC数据集确定火山电离层TEC异常扰动的上界线、下界线。

[0170] 在一些实施例中，通过Prophet算法改进LSTM算法，基于最优电离层TEC数据集确定火山电离层TEC异常扰动的上界线、下界线，具体包括：

[0171] 分别赋予Prophet算法的预测值 和LSTM模型预测值 的i时刻电离层TEC的初始权重ω1、ω2。

[0172] 通过与实测电离层TEC进行误差精度分析，确定目标函数：

[0173]

[0174] 其中，N表示探测时段的总数；yi表示实测电离层TEC； 表示i时刻Prophet算法的TEC预测值， 表示i时刻LSTM算法的TEC预测值。

[0175] 预设约束条件：ω1+ω2＝1,ω1≥0,ω2≥0；

[0176] 利用最小二乘法确定最优权重ω1和ω2，并确定火山电离层TEC异常扰动的上界线、下界线。

[0177] 具体的，公式(14)中的y(t)即Prophet算法的预测值 公式(15)中的ht即LSTM模型预测值

[0178] 在一些实施例中，利用最小二乘法确定最优权重ω1和ω2，并确定火山电离层TEC异常扰动的上界线、下界线，具体包括：

[0179] Error＝TECactual‑TECpredicted  (17)

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

[0181] LB＝Mean(Error)‑2×Std(Error)  (19)

[0182] 式中，TECactual表示实测电离层TEC；TECpredicted表示新算法预测的电离层TEC；Mean表示取平均值；Std表示求标准差；UB表示电离层TEC异常的上界线；LB表示电离层TEC异常的下界线。

[0183] S105、根据电离层TEC观测数据与火山电离层TEC异常扰动的上界线、下界线的比较关系，确定火山电离层是否发生异常扰动。

[0184] 通过将电离层TEC观测数据与确定的火山电离层TEC异常扰动的上界线、下界线进行比较，可以清晰地识别出火山活动是否引起了电离层的异常扰动，这种方法不仅科学且高效，能够为火山监测和预警提供重要参考。

[0185] 在一些实施例中，ω1和ω2的值分别设为0.6和0.4，图4和图5为利用本方法探测的2022年汤加火山前电离层异常情况，如图4所示为：2022年1月1日—16日SAMO站上空的电离层TEC变化(红色线表示TEC真实值，蓝色线表示算法预测值，红色粗实线表示火山喷发时刻)，红色粗实线标记了火山喷发的具体时刻，通过对比红蓝两线，可以清晰地看到在火山喷发前后电离层TEC的显著波动。在火山喷发前，算法预测值(蓝色线)与真实值(红色线)保持较为接近的趋势，但在火山喷发时刻及其后的一段时间内，预测值出现了较大的偏离，反映了火山喷发对电离层的强烈扰动。

[0186] 如图5所示为：2022年1月1日—16日SAMO站上空的电离层TEC异常(黑色线和蓝色线分别表示TEC异常的上下界线，红色粗实线表示火山喷发时刻)；红色粗实线同样标记了火山喷发的时刻，从图中可以看到，在火山喷发前后，TEC值迅速跨越了异常界限，表明火山活动对电离层产生了显著的异常扰动，通过对比TEC值与异常界限的偏离程度，可以量化评估火山喷发对电离层造成的具体影响程度和持续时间。

[0187] 综上，在火山喷发前，本申请提出的方法能够较好地预测电离层TEC的变化趋势，通过设定TEC异常的上下界限，有效地识别出了火山喷发引起的电离层异常扰动，证明了该方案在火山电离层异常扰动的确定上具有较高的灵敏度和准确性，通过实时监测电离层TEC的变化，可以提前发现火山活动的迹象，为相关部门提供宝贵的预警时间。

[0188] 本申请还提出了一种火山电离层异常扰动探测系统，系统包括：数据获取单元、判断条件预设单元和确定单元；

[0189] 数据获取单元，用于获取地基GNSS原始观测数据；对地基GNSS原始观测数据进行预处理，获取电离层TEC数据集；根据电离层TEC数据集的小波能量谱测度时间序列在时域和频域内的方差，去除与空间天气环境相关性高的电离层TEC数据，获取不受空间天气环境影响的最优电离层TEC数据集；

[0190] 判断条件预设单元，用于通过Prophet算法改进LSTM算法，基于最优电离层TEC数据集确定火山电离层TEC异常扰动的上界线、下界线；

[0191] 确定单元，用于根据电离层TEC数据集与火山电离层TEC异常扰动的上界线、下界线的比较关系，确定火山电离层是否发生异常扰动。

[0192] 本申请还提出了一种可读存储介质，存储有计算机程序，计算机程序被处理器执行时，使得处理器执行如下步骤：

[0193] 获取地基GNSS原始观测数据；

[0194] 对地基GNSS原始观测数据进行预处理，获取电离层TEC数据集；

[0195] 根据电离层TEC数据集的小波能量谱测度时间序列在时域和频域内的方差，去除与空间天气环境相关性高的电离层TEC数据，获取不受空间天气环境影响的最优电离层TEC数据集；

[0196] 通过Prophet算法改进LSTM算法，基于最优电离层TEC数据集确定火山电离层TEC异常扰动的上界线、下界线；

[0197] 根据电离层TEC观测数据与火山电离层TEC异常扰动的上界线、下界线的比较关系，确定火山电离层是否发生异常扰动。

[0198] 本申请还提出了一种计算机设备，包括存储器和处理器，存储器存储有计算机程序，计算机程序被处理器执行如下步骤：

[0199] 获取地基GNSS原始观测数据；

[0200] 对地基GNSS原始观测数据进行预处理，获取电离层TEC数据集；

[0201] 根据电离层TEC数据集的小波能量谱测度时间序列在时域和频域内的方差，去除与空间天气环境相关性高的电离层TEC数据，获取不受空间天气环境影响的最优电离层TEC数据集；

[0202] 通过Prophet算法改进LSTM算法，基于最优电离层TEC数据集确定火山电离层TEC异常扰动的上界线、下界线；

[0203] 根据电离层TEC观测数据与火山电离层TEC异常扰动的上界线、下界线的比较关系，确定火山电离层是否发生异常扰动。

[0204] 本领域普通技术人员能够理解实现上述实施例方法中的全部或部分流程，是能够通过计算机程序来指令相关的硬件来完成，所述的程序可存储于一非易失性计算机可读取存储介质中，该程序在执行时，可包括如上述各方法的实施例的流程。其中，本申请所提供的各实施例中所使用的对存储器、存储、数据库或其它介质的任何引用，均可包括非易失性和/或易失性存储器。非易失性存储器可包括只读存储器(ROM)、可编程ROM(PROM)、电可编程ROM(EPROM)、电可擦除可编程ROM(EEPROM)或闪存。易失性存储器可包括随机存取存储器(RAM)或者外部高速缓冲存储器。作为说明而非局限，RAM以多种形式可得，诸如静态RAM(SRAM)、动态RAM(DRAM)、同步DRAM(SDRAM)、双数据率SDRAM(DDRSDRAM)、增强型SDRAM(ESDRAM)、同步链路(Synchlink)DRAM(SLDRAM)、存储器总线(Rambus)直接RAM(RDRAM)、直接存储器总线动态RAM(DRDRAM)、以及存储器总线动态RAM(RDRAM)等。

[0205] 以上实施例的各技术特征能够进行任意的组合，为使描述简洁，未对上述实施例中的各个技术特征所有可能的组合都进行描述，然而，只要这些技术特征的组合不存在矛盾，都应当认为是本说明书记载的范围。

[0206] 以上所述实施例仅表达了本申请的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对本申请专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本申请构思的前提下，还能够做出若干变形和改进，这些都属于本申请的保护范围，以上所揭露的仅为本发明较佳实施例而已，当然不能以此来限定本发明之权利范围，因此依本发明权利要求所作的等同变化，仍属本发明所涵盖的范围。