CN116299595A - 面向扰动探测的残差电离层层析方法、装置及介质 - Google Patents
面向扰动探测的残差电离层层析方法、装置及介质 Download PDFInfo
- Publication number
- CN116299595A CN116299595A CN202310517850.8A CN202310517850A CN116299595A CN 116299595 A CN116299595 A CN 116299595A CN 202310517850 A CN202310517850 A CN 202310517850A CN 116299595 A CN116299595 A CN 116299595A
- Authority
- CN
- China
- Prior art keywords
- value
- stec
- ionosphere
- fitting
- vtec
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 239000005433 ionosphere Substances 0.000 title claims abstract description 95
- 238000000034 method Methods 0.000 title claims abstract description 83
- 238000004587 chromatography analysis Methods 0.000 title claims abstract description 27
- 238000001514 detection method Methods 0.000 title claims description 18
- 230000008859 change Effects 0.000 claims abstract description 33
- 230000002159 abnormal effect Effects 0.000 claims abstract description 18
- 238000011160 research Methods 0.000 claims description 25
- 230000008569 process Effects 0.000 claims description 14
- 230000007812 deficiency Effects 0.000 claims description 11
- 238000005516 engineering process Methods 0.000 claims description 9
- 238000012937 correction Methods 0.000 claims description 8
- 238000005259 measurement Methods 0.000 claims description 8
- 238000009499 grossing Methods 0.000 claims description 7
- 239000011159 matrix material Substances 0.000 claims description 4
- 230000006870 function Effects 0.000 description 12
- 238000010586 diagram Methods 0.000 description 6
- 230000000694 effects Effects 0.000 description 4
- 238000011084 recovery Methods 0.000 description 4
- 230000005856 abnormality Effects 0.000 description 3
- 230000007717 exclusion Effects 0.000 description 3
- 238000003325 tomography Methods 0.000 description 3
- 208000028257 Joubert syndrome with oculorenal defect Diseases 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 2
- 230000007547 defect Effects 0.000 description 2
- 230000006978 adaptation Effects 0.000 description 1
- 230000004075 alteration Effects 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 238000004364 calculation method Methods 0.000 description 1
- 238000004891 communication Methods 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 239000003574 free electron Substances 0.000 description 1
- 230000005358 geomagnetic field Effects 0.000 description 1
- 230000010354 integration Effects 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 230000000116 mitigating effect Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000005855 radiation Effects 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 230000003313 weakening effect Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
- G01S19/01—Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
- G01S19/13—Receivers
- G01S19/35—Constructional details or hardware or software details of the signal processing chain
- G01S19/37—Hardware or software details of the signal processing chain
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
- G01S19/01—Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
- G01S19/03—Cooperating elements; Interaction or communication between different cooperating elements or between cooperating elements and receivers
- G01S19/07—Cooperating elements; Interaction or communication between different cooperating elements or between cooperating elements and receivers providing data for correcting measured positioning data, e.g. DGPS [differential GPS] or ionosphere corrections
- G01S19/072—Ionosphere corrections
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
- G01S19/01—Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
- G01S19/13—Receivers
- G01S19/23—Testing, monitoring, correcting or calibrating of receiver elements
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
- Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation
Landscapes
- Engineering & Computer Science (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Computer Networks & Wireless Communication (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Signal Processing (AREA)
- Position Fixing By Use Of Radio Waves (AREA)
Abstract
本发明属于GNSS电离层反演领域电离层层析技术领域，实现反演电离层扰动发生时电子密度变化的三维时空分布，解决电子密度异常变化量的三维反演、地面观测站分布不均、迭代算法未结合实际电子分布情况等问题。具体公开一种面向扰动探测的残差电离层层析方法、装置及介质，采用布设虚拟观测站得到虚拟观测值来解决，在反演区域站点稀疏的地区布设虚拟观测站，根据虚拟站与GNSS卫星连线计算虚拟穿刺点，利用实际观测数据建立区域电离层模型或GIM内插方法得到虚拟穿刺点的虚拟观测值。
Description
技术领域
本发明属于GNSS电离层反演领域电离层层析技术领域，更具体地，涉及一种面向扰动探测的残差电离层层析方法、装置及介质。
背景技术
电离层是近地空间环境中直接影响人类生活的重要一环，其自身处于复杂多变的日地空间环境中，经常发生结构偏离常规形态的扰动，扰动发生时会对无线电通信、卫星导航定位、航天器安全等产生严重影响。随着全球导航卫星系统(Global NavigationSatellite System，GNSS)逐渐发展完善，其在电离层研究领域的作用逐渐发掘。GNSS电离层层析技术基于导航卫星电磁波信号延迟来反演电离层电子密度空间分布，经过近三十年的发展与完善已经日臻成熟，凭借其全天时全天候不间断观测、能够得到高精度高时空分辨率的三维电离层电子密度信息等优点在电离层扰动研究中应用广泛。
电离层层析是典型的秩亏问题，因为存在如GNSS观测站地面分布不均匀，划定反演区域内观测站过于稀少，GNSS卫星观测视角有限等诸多问题，同时层析射线还要求必须同时穿过电离层格网的顶层和底层以保证STEC的准确性，种种限制均会导致许多格网点没有射线穿过，得不到任何观测信息。层析反演的过程就是求层析方程组的解得到格网内的电子密度。目前像素基层析中常用的迭代重构算法有代数重构算法、乘法代数重构算法、同时迭代重构算法等。
不过传统层析方法得到的结果是反演时刻电离层的电子密度，无法判断电离层是否发生扰动。利用电离层总电子含量(Total Electron Content，TEC))相对背景场的残差建立层析方程，能够实现反演电离层电子密度变化的时空分布。
Zhai等[Zhai C,Yao Y,Kong J.Three-dimensional reconstruction ofseismo-traveling ionospheric disturbances after March 11,2011,Japan Tohokuearthquake[J].Journal of Geodesy,2021,95(7).]基于穿刺点轨迹求电离层残差，假设卫星穿刺点轨迹每天相同，穿刺点的倾斜路径总电子含量(Slant Total ElectronContent,STEC)在电离层正常变化情况下是平滑的，取当前时刻前后5分钟作为窗口，利用滑动窗口取平均值去趋势求STEC残差，建立层析方程，利用同时代数重构技术计算层析结果。该方法可以有效探测出短期内电离层的突然扰动。不过探测效果与滑动窗口的选取有关，窗口过小无法探测出异常，窗口过大时背景STEC变化趋势与穿刺点处STEC物理关联性减弱；直接采用平均值的方式计算背景值忽略了STEC自身的变化规律；同时当电离层扰动持续时间大于探测窗口乃至达到数小时时，其扰动本身也会被当做趋势项去除，本方法不适用于长时间、大范围的电离层扰动研究。
He[He L,Heki K.Three-dimensional distribution of ionosphericanomalies prior to three large earthquakes in Chile[J].Geophysical ResearchLetters,2016,43(14):7287-93]，[He L,Heki K.Three-Dimensional Tomography ofIonospheric Anomalies Immediately Before the 2015Illapel Earthquake,CentralChile[J].Journal of Geophysical Research:Space Physics,2018,123(5):4015-25.]与Heki[Heki K.Ionospheric electron enhancement preceding the 2011Tohoku-Okiearthquake[J].Geophysical Research Letters,2011,38(17):n/a-n/a.]、[Heki K,Enomoto Y.Mw dependence of the preseismic ionospheric electron enhancements[J].Journal of Geophysical Research:Space Physics,2015,120(8):7006-20.]等利用参考曲线求电离层残差，基于电离层天顶方向总电子含量(Vertical Total ElectronContent,VTEC)的时间序列用多项式拟合建立以时间为变量的参考曲线，通过实际观测值与参考曲线比较作差，得到VTEC异常变化，再由投影函数转换成STEC的异常变化，完成层析计算。在地震电离层扰动研究中，为避免地震前后对电离层的影响，设立了从震前40分钟至震后20分钟的“排除时间”。本方法时间序列的开始时间(去掉“排除时间”后的开始时间)没有得到客观约束；而且建立参考曲线时引入了震后观测数据，虽然“排除时间”包含震后20分钟，但仍然难以保证震后电离层异常不会对参考曲线产生偏移；同时也正是“排除时间”的存在决定了该方法不适用于临震短期电离层异常探测。此外，该方法同样不适用于持续时间较长的电离层扰动研究。
电离层层析中GNSS信号射线在像素格网中的分布高度依赖于GNSS观测站的布设，可是GNSS观测站点的实际分布情况并不理想，在中国、北美、欧洲等地区站点分布密集，而在非洲、中东等地区以及沙漠、海洋地区站点稀疏甚至没有观测站。严重影响了层析技术的应用范围和计算结果精度。
发明内容
提供了本发明以解决现有技术中存在的上述问题。因此，需要一种面向扰动探测的残差电离层层析方法、装置及介质，实现反演电离层扰动发生时电子密度变化的三维时空分布，解决电子密度异常变化量的三维反演、地面观测站分布不均、迭代算法未结合实际电子分布情况等问题。
根据本发明的第一方案，提供了一种面向扰动探测的残差电离层层析方法，所述方法包括：
获取观测站数据，所述观测站数据包括观测站的观测值文件、GPS导航电文文件、精密星历和差分码偏差文件，根据地面站的GPS观测数据计算STEC的观测值；
基于穿刺点时间序列求电离层残差，采用区域电离层建模或GIM数据内插方法计算该穿刺点在同一时刻不同日期的VTEC，得到穿刺点VTEC时间序列后，通过多项式拟合计算研究时刻的拟合值，并通过投影函数得到STEC的拟合值，将STEC的拟合值与实测值比较得到STEC的异常变化；
在层析区域建立格网模型，令每个体格网中的电子密度相同，建立残差电离层层析方程；
利用同时代数重构技术解决层析方程组秩亏问题求解未知量。
进一步地，所述根据地面站的GPS观测数据计算STEC的观测值，具体包括：
通过GPS的双频信号计算得到所述STEC的观测值：
式中，f1和f2是GPS信号的频率，f1＝1575.42MHz,f2＝1227.60MHz,
和/>
是f1和f2频率信号的载波相位平滑伪距值；Δbk和Δbs分别是接收机和卫星的差分码偏差。
进一步地，基于穿刺点时间序列求电离层残差，采用区域电离层建模或GIM数据内插方法计算该穿刺点在同一时刻不同日期的VTEC，得到穿刺点VTEC时间序列后，通过多项式拟合计算研究时刻的拟合值，并通过投影函数得到STEC的拟合值，将STEC的拟合值与实测值比较得到STEC的异常变化的具体表达式为：
式中，
为通过区域建模或GIM内插得到的ti时刻穿刺点VTEC值；VTEC(β,s)和VTEC(β,λ)分别是电离层区域模型和GIM内插，β为地理纬度，s为太阳时角，λ为地理经度，
为研究时刻t0的穿刺点VTEC拟合值；an为待求参数；/>
是穿刺点处STEC实测值；f(Z)是投影函数；/>
为求得穿刺点的STEC残差；Z为测站处卫星天顶角，Z′为穿刺点处卫星天顶角。
进一步地，根据研究范围确定采用区域电离层建模或GIM数据内插方法。
进一步地，所述残差电离层层析方程表示为：
dSTECm×1＝Am×n·Δxn×1+εm×1 (3)
式中，dSTEC是GNSS卫星信号传播路径上电离层STEC残差组成的m维列向量，m是测量次数，n是反演区域的体素格网个数，各条射线在对应格网中的截距组成m×n阶矩阵A，Δx为格网内电子密度的变化量，ε为随机噪声。
进一步地，利用同时代数重构技术解决层析方程组秩亏问题求解未知量的具体表达式如下：
式中，
是第j个格网经过k+1次迭代后的电子密度变化值；P是第j个格网内经过的射线数；λ是松弛因子，用于控制收敛性，取值为0<λ<1；yi是第i条射线的dSTEC值；Ai,j是第i条射线在第j个格网中的截距；IEDj是第j个格网的电子密度IRI模型值；D是第i条射线上的修正量；W是第i条射线上第j个格网的dSTEC修正权重。
根据本发明的第二技术方案，提供一种面向扰动探测的残差电离层层析方法装置，所述装置包括：
获取单元，被配置为获取观测站数据，所述观测站数据包括观测站的观测值文件、GPS导航电文文件、精密星历和差分码偏差文件，根据地面站的GPS观测数据计算STEC的观测值；
比较单元，被配置为基于穿刺点时间序列求电离层残差，采用区域电离层建模或GIM数据内插方法计算该穿刺点在同一时刻不同日期的VTEC，得到穿刺点VTEC时间序列后，通过多项式拟合计算研究时刻的拟合值，并通过投影函数得到STEC的拟合值，将STEC的拟合值与实测值比较得到STEC的异常变化；
建立单元，被配置为在层析区域建立格网模型，令每个体格网中的电子密度相同，建立残差电离层层析方程；
求解单元，被配置为利用同时代数重构技术解决层析方程组秩亏问题求解未知量。
进一步地，所述获取单元被进一步配置为通过GPS的双频信号计算得到所述STEC的观测值：
式中，f1和f2是GPS信号的频率，f1＝1575.42MHz,f2＝1227.60MHz,
和/>
是f1和f2频率信号的载波相位平滑伪距值；Δbk和Δbs分别是接收机和卫星的差分码偏差。
进一步地，所述比较单元被进一步配置为基于穿刺点时间序列求电离层残差，采用区域电离层建模或GIM数据内插方法计算该穿刺点在同一时刻不同日期的VTEC，得到穿刺点VTEC时间序列后，通过多项式拟合计算研究时刻的拟合值，并通过投影函数得到STEC的拟合值，将STEC的拟合值与实测值比较得到STEC的异常变化的具体表达式为：
式中，
为通过区域建模或GIM内插得到的ti时刻穿刺点VTEC值；VTEC(β,s)和VTEC(β,λ)分别是电离层区域模型和GIM内插，β为地理纬度，s为太阳时角，λ为地理经度，
为研究时刻t0的穿刺点VTEC拟合值；an为待求参数；/>
是穿刺点处STEC实测值；f(Z)是投影函数；/>
为求得穿刺点的STEC残差；Z为测站处卫星天顶角，Z′为穿刺点处卫星天顶角。
根据本发明的第三方面，提供一种可读存储介质，所述可读存储介质存储有一个或者多个程序，所述一个或者多个程序可被一个或者多个处理器执行，以实现如上所述的方法。
本发明至少具有以下有益效果：
本发明实现探测分析长时间、大范围的电离层扰动特征，扰动本身不会被当作变化趋势被去除，克服站点分布不均匀的缺陷，适用于磁暴、地震等多种因素导致的电离层扰动。为电离层扰动研究提供大量的起止、传播、变化等信息，为减轻电离层扰动的危害、探索地震电离层异常机理做出贡献。
附图说明
图1为本发明实施例的一种面向扰动探测的残差电离层层析方法的完整运算流程图；
图2为本发明实施例的一种面向扰动探测的残差电离层层析方法的流程图；
图3为本发明实施例的一种面向扰动探测的残差电离层层析装置的结构图；
图4为本发明实施例的选用GNSS观测站全球分布图；
图5为本发明实施例的磁暴期间电离层扰动电子密度三维变化情况示意图；
图6为本发明实施例的磁暴期间375km高度电离层电子密度变化情况示意图；
图7为本发明实施例的磁暴期间电离层全高度积分TEC变化情况示意图；
图8为本发明实施例的磁暴期间东经100.5°电离层电子密度变化剖面图；
图9为本发明实施例的磁暴期间西经136.5°电离层电子密度变化剖面图。
具体实施方式
为使本领域技术人员更好的理解本发明的技术方案，下面结合附图和具体实施方式对本发明作详细说明。下面结合附图和具体实施例对本发明的实施例作进一步详细描述，但不作为对本发明的限定。本文中所描述的各个步骤，如果彼此之间没有前后关系的必要性，则本文中作为示例对其进行描述的次序不应视为限制，本领域技术人员应知道可以对其进行顺序调整，只要不破坏其彼此之间的逻辑性导致整个流程无法实现即可。
本发明实施例提供一种面向扰动探测的残差电离层层析方法，请参阅图1所示，该方法在层析的实际运算过程中，首先利用虚拟观测值按照基于穿刺点时间序列求得虚拟电离层残差，进行残差电离层层析；而后在计算结果的基础上利用实际观测值按同样方法进行残差电离层层析。
具体来说，请参阅图2，该方法包括如下步骤：
步骤S100，获取观测站数据，所述观测站数据包括观测站的观测值文件、GPS导航电文文件、精密星历和差分码偏差文件，根据地面站的GPS观测数据计算STEC的观测值。
需要说明的是，观测站数据从各地的连续运行参考站网(ContinuouslyOperating Reference Stations，CORS)和国际GNSS服务中心(International GNSSService,IGS)获取。
在一些实施例中，所述根据地面站的GPS观测数据计算STEC的观测值，具体包括：
通过GPS的双频信号计算得到所述STEC的观测值：
式中，f1和f2是GPS信号的频率，f1＝1575.42MHz,f2＝1227.60MHz,
和/>
是f1和f2频率信号的载波相位平滑伪距值；Δbk和Δbs分别是接收机和卫星的差分码偏差。
步骤S200，基于穿刺点时间序列求电离层残差，采用区域电离层建模或GIM数据内插方法计算该穿刺点在同一时刻不同日期的VTEC，得到穿刺点VTEC时间序列后，通过多项式拟合计算研究时刻的拟合值，并通过投影函数得到STEC的拟合值，将STEC的拟合值与实测值比较得到STEC的异常变化。
在一些实施例中，根据研究范围确定采用区域电离层建模或GIM数据内插方法。当研究范围较大例如全球时，采用GIM内插方法，当研究区域较小且精度要求较高时利用实测数据建立电离层区域模型，建模或内插的时间范围为研究时刻前地磁活动平静的16天。
在一些实施例中，基于穿刺点时间序列求电离层残差，采用区域电离层建模或GIM数据内插方法计算该穿刺点在同一时刻不同日期的VTEC，得到穿刺点VTEC时间序列后，通过多项式拟合计算研究时刻的拟合值，并通过投影函数得到STEC的拟合值，将STEC的拟合值与实测值比较得到STEC的异常变化的具体表达式为：
式中，
为通过区域建模或GIM内插得到的ti时刻穿刺点VTEC值；VTEC(β,s)和VTEC(β,λ)分别是电离层区域模型和GIM内插，β为地理纬度，s为太阳时角，λ为地理经度，
为研究时刻t0的穿刺点VTEC拟合值；an为待求参数；/>
是穿刺点处STEC实测值；f(Z)是投影函数；/>
为求得穿刺点的STEC残差；Z为测站处卫星天顶角，Z′为穿刺点处卫星天顶角。
步骤S300，在层析区域建立格网模型，令每个体格网中的电子密度相同，建立残差电离层层析方程。
在一些实施例中，所述残差电离层层析方程表示为：
dSTECm×1＝Am×n·Δxn×1+εm×1 (3)
式中，dSTEC是GNSS卫星信号传播路径上电离层STEC残差组成的m维列向量，m是测量次数，n是反演区域的体素格网个数，各条射线在对应格网中的截距组成m×n阶矩阵A，Δx为格网内电子密度的变化量，ε为随机噪声。
步骤S400，利用同时代数重构技术解决层析方程组秩亏问题求解未知量。
在一些实施例中，利用同时代数重构技术解决层析方程组秩亏问题求解未知量的具体表达式如下：
式中，
是第j个格网经过k+1次迭代后的电子密度变化值；P是第j个格网内经过的射线数；λ是松弛因子，用于控制收敛性，取值为0<λ<1；yi是第i条射线的dSTEC值；Ai,j是第i条射线在第j个格网中的截距；IEDj是第j个格网的电子密度IRI模型值；D是第i条射线上的修正量；W是第i条射线上第j个格网的dSTEC修正权重。
本发明实施例还提供一种面向扰动探测的残差电离层层析方法装置，如图3所示，所述装置300包括：
获取单元301，被配置为获取观测站数据，所述观测站数据包括观测站的观测值文件、GPS导航电文文件、精密星历和差分码偏差文件，根据地面站的GPS观测数据计算STEC的观测值；
比较单元302，被配置为基于穿刺点时间序列求电离层残差，采用区域电离层建模或GIM数据内插方法计算该穿刺点在同一时刻不同日期的VTEC，得到穿刺点VTEC时间序列后，通过多项式拟合计算研究时刻的拟合值，并通过投影函数得到STEC的拟合值，将STEC的拟合值与实测值比较得到STEC的异常变化；
建立单元303，被配置为在层析区域建立格网模型，令每个体格网中的电子密度相同，建立残差电离层层析方程；
求解单元304，被配置为利用同时代数重构技术解决层析方程组秩亏问题求解未知量。
在一些实施例中，所述获取单元被进一步配置为通过GPS的双频信号计算得到所述STEC的观测值：
式中，f1和f2是GPS信号的频率，f1＝1575.42MHz,f2＝1227.60MHz,
和/>
是f1和f2频率信号的载波相位平滑伪距值；Δbk和Δbs分别是接收机和卫星的差分码偏差。
在一些实施例中，所述比较单元被进一步配置为基于穿刺点时间序列求电离层残差，采用区域电离层建模或GIM数据内插方法计算该穿刺点在同一时刻不同日期的VTEC，得到穿刺点VTEC时间序列后，通过多项式拟合计算研究时刻的拟合值，并通过投影函数得到STEC的拟合值，将STEC的拟合值与实测值比较得到STEC的异常变化的具体表达式为：
式中，
为通过区域建模或GIM内插得到的ti时刻穿刺点VTEC值；VTEC(β,s)和VTEC(β,λ)分别是电离层区域模型和GIM内插，β为地理纬度，s为太阳时角，λ为地理经度，
为研究时刻t0的穿刺点VTEC拟合值；an为待求参数；/>
是穿刺点处STEC实测值；f(Z)是投影函数；/>
为求得穿刺点的STEC残差；Z为测站处卫星天顶角，Z′为穿刺点处卫星天顶角。
在一些实施例中，所述比较单元被进一步配置为根据研究范围确定采用区域电离层建模或GIM数据内插方法。
在一些实施例中，所述残差电离层层析方程表示为：
dSTECm×1＝Am×n·Δxn×1+εm×1 (3)
式中，dSTEC是GNSS卫星信号传播路径上电离层STEC残差组成的m维列向量，m是测量次数，n是反演区域的体素格网个数，各条射线在对应格网中的截距组成m×n阶矩阵A，Δx为格网内电子密度的变化量，ε为随机噪声。
在一些实施例中，所述求解单元被进一步配置为利用同时代数重构技术解决层析方程组秩亏问题求解未知量的具体表达式如下：
式中，
是第j个格网经过k+1次迭代后的电子密度变化值；P是第j个格网内经过的射线数；λ是松弛因子，用于控制收敛性，取值为0<λ<1；yi是第i条射线的dSTEC值；Ai,j是第i条射线在第j个格网中的截距；IEDj是第j个格网的电子密度IRI模型值；D是第i条射线上的修正量；W是第i条射线上第j个格网的dSTEC修正权重。
需要说明的是，本实施例所提出的装置与在先阐述的方法属于同一技术思路，其局有限相同的工作原理，并能起到同样的有益效果，此处不赘述。
本发明实施例提供一种可读存储介质，所述可读存储介质存储有一个或者多个程序，所述一个或者多个程序可被一个或者多个处理器执行，以实现如上任一实施例所述的面向扰动探测的残差电离层层析方法。
利用本申请的技术方案，基于全球7个CORS网2143个观测站的观测数据，实现2018年8月25-26日磁暴期间全球尺度电离层扰动的层析反演，技术效果如下。
图4为选用GNSS观测站全球分布图。
本实施例划分格网分辨率为经度3°，纬度1.5°，高程50km，研究范围为纬度87°N～87°S，经度180°W～180°E，高程100～1000km，时间分辨率为1小时，卫星截止高度角15°。在全球范围布设纬度1.5°，经度3°分辨率的虚拟观测站，根据精密星历文件提供GNSS卫星位置信息，鉴于研究尺度为全球范围同时研究对象为磁暴期间电离层扰动，采用GIM数据内插方法得到各虚拟穿刺点的VTEC值，以弥补观测数据的缺失。
图5展示磁暴期间电离层扰动的层析结果，分别为高度剖面和经度剖面，x轴表示地理经度，y轴表示纬度，z轴表示高程。可以发现电子密度变化主要发生在高度100km-600km的电离层中且集中于300km-450km处，600km以上高度的电子变化不明显，这与平静时期电离层中电子密度的实际分布情况有关。在地磁活动趋于最剧烈的8月26日06:00UTC，电子密度变化也达到最大，在100km高程和600km高程以上也可以发现电子密度的增加。磁暴期间电离层扰动的范围遍及全球，不过赤道和低纬地区更加明显，并且存在赤道电离层异常现象，电子密度变化量远大于3×105el/cm3。在磁暴恢复相期间，电子密度变化在南北半球呈现明显的不对称现象，北半球电子密度降低而南半球电子密度变化值依然增加
图6介绍磁暴期间2018年8月25日17:00UTC至8月26日20:00UTC375km高度电离层电子密度变化情况，红线为磁赤道线。可以发现8月25日17:00UTC已经有电子密度增加的趋势，磁暴初始相阶段北极和欧洲北部地区存在电子密度降低的现象，电子密度变化量可达-0.87～-0.11×105el/cm3，至磁暴正相阶段电离层扰动范围已经遍及全球，电子密度变化量达到1.47～8.55×105el/cm3。同时在磁赤道南北两侧存在赤道双峰现象。在磁暴恢复相阶段电离层电子密度变化情况呈现南北不对称现象，北半球电离层呈现负异常而南半球电离层依然呈现正异常，这种现象一直持续到电离层恢复平静。虽然电离层扰动范围随时间变化有所不同，不过整体趋势表现为自东向西移动，这应该与地球自转有关，地球自转导致太阳辐射作用在地球的位置发生变化。
图7为电离层电子密度变化量在各个高度积分后的结果，可以发现其变化趋势和电离层扰动特征基本和375km高程处电子密度变化情况保持一致。注意到在磁暴恢复相阶段电离层扰动南北半球不对称时，北半球电子密度降低不论刚刚显现还是非常显著，其扰动分布范围与南半球的电子密度增加分布范围在磁赤道附近有明显的分界。因此推测导致南北半球不对称的原因除了热层夏冬季节环流外，磁暴恢复时地磁场的迅速变化也对自由电子的传播运动产生重要影响。
图8表现磁暴期间电离层电子密度变化情况在东经100.5°的剖面图，可以发现本地区电离层扰动从2018年8月25日23:00UTC开始逐渐剧烈，扰动范围在纬度向从北纬24N～50N、南纬0～30S逐步扩散至中低纬大部分地区，高度范围从200km～400km扩散至100km～600km，电子密度变化值从0.68～1.79×105el/cm3增加到1.59～7.75×105el/cm3，从在赤道南北两侧存在明显的双峰现象。且在北极存在负异常，与前文结论相符。在强烈扰动在持续约4小时后，从8月26日11:00UTC开始电离层扰动开始减弱，在扰动平复过程中电子密度增加值逐渐较少，之后在大部分纬向范围电子密度变化值为负，南北半球不对称现象虽然存在但是并不明显。
图9展示磁暴期间电离层电子密度变化情况在西经136.5°的剖面图，可以看出在2018年8月25日17:00UTC时本区域电离层扰动已经非常明显，至8月25日23:00UTC已经达到峰值，扰动纬度向范围从北纬30N～47N和52S～6N扩散至全纬度，高度范围也从150km～350km扩散至100km～600km，电子密度变化值从0.39～1.19×105el/cm3增加到1.68～4.09×105el/cm3，在磁赤道南北两侧存在明显的双峰现象。在电离层扰动减弱期间，南北半球均出现电离层负异常现象，8月26日15:00UTC后出现明显的南北半球不对称现象。
此外，尽管已经在本文中描述了示例性实施例，其范围包括任何和所有基于本发明的具有等同元件、修改、省略、组合(例如，各种实施例交叉的方案)、改编或改变的实施例。权利要求书中的元件将被基于权利要求中采用的语言宽泛地解释，并不限于在本说明书中或本申请的实施期间所描述的示例，其示例将被解释为非排他性的。因此，本说明书和示例旨在仅被认为是示例，真正的范围和精神由以下权利要求以及其等同物的全部范围所指示。
以上描述旨在是说明性的而不是限制性的。例如，上述示例(或其一个或更多方案)可以彼此组合使用。例如本领域普通技术人员在阅读上述描述时可以使用其它实施例。另外，在上述具体实施方式中，各种特征可以被分组在一起以简单化本发明。这不应解释为一种不要求保护的发明的特征对于任一权利要求是必要的意图。相反，本发明的主题可以少于特定的发明的实施例的全部特征。从而，以下权利要求书作为示例或实施例在此并入具体实施方式中，其中每个权利要求独立地作为单独的实施例，并且考虑这些实施例可以以各种组合或排列彼此组合。本发明的范围应参照所附权利要求以及这些权利要求赋权的等同形式的全部范围来确定。
Claims (10)
1.一种面向扰动探测的残差电离层层析方法，其特征在于，所述方法包括：
获取观测站数据，所述观测站数据包括观测站的观测值文件、GPS导航电文文件、精密星历和差分码偏差文件，根据地面站的GPS观测数据计算STEC的观测值；
基于穿刺点时间序列求电离层残差，采用区域电离层建模或GIM数据内插方法计算该穿刺点在同一时刻不同日期的VTEC，得到穿刺点VTEC时间序列后，通过多项式拟合计算研究时刻的拟合值，并通过投影函数得到STEC的拟合值，将STEC的拟合值与实测值比较得到STEC的异常变化；
在层析区域建立格网模型，令每个体格网中的电子密度相同，建立残差电离层层析方程；
利用同时代数重构技术解决层析方程组秩亏问题求解未知量。方法整体分为两步，首先利用虚拟观测值按照基于穿刺点时间序列求得虚拟电离层残差，进行残差电离层层析；而后在计算结果的基础上利用实际观测值按同样方法进行残差电离层层析。
2.根据权利要求1所述的方法，其特征在于，所述根据地面站的GPS观测数据计算STEC的观测值，具体包括：
通过GPS的双频信号计算得到所述STEC的观测值：
式中，f1和f2是GPS信号的频率，f1＝1575.42MHz,f2＝1227.60MHz,
和/>
是f1和f2频率信号的载波相位平滑伪距值；Δbk和Δbs分别是接收机和卫星的差分码偏差。
3.根据权利要求1所述的方法，其特征在于，基于穿刺点时间序列求电离层残差，采用区域电离层建模或GIM数据内插方法计算该穿刺点在同一时刻不同日期的VTEC，得到穿刺点VTEC时间序列后，通过多项式拟合计算研究时刻的拟合值，并通过投影函数得到STEC的拟合值，将STEC的拟合值与实测值比较得到STEC的异常变化的具体表达式为：
式中，
为通过区域建模或GIM内插得到的ti时刻穿刺点VTEC值；VTEC(β,s)和VTEC(β,λ)分别是电离层区域模型和GIM内插，β为地理纬度，s为太阳时角，λ为地理经度，
为研究时刻t0的穿刺点VTEC拟合值；an为待求参数；/>
是穿刺点处STEC实测值；f(Z)是投影函数；/>
为求得穿刺点的STEC残差；Z为测站处卫星天顶角，Z′为穿刺点处卫星天顶角。
4.根据权利要求3所述的方法，其特征在于，根据研究范围确定采用区域电离层建模或GIM数据内插方法。
5.根据权利要求1所述的方法，其特征在于，所述残差电离层层析方程表示为：
dSTECm×1＝Am×n.Δxn×1+εm×1 (3)
式中，dSTEC是GNSS卫星信号传播路径上电离层STEC残差组成的m维列向量，m是测量次数，n是反演区域的体素格网个数，各条射线在对应格网中的截距组成m×n阶矩阵A，Δx为格网内电子密度的变化量，ε为随机噪声。
6.根据权利要求1所述的方法，其特征在于，利用同时代数重构技术解决层析方程组秩亏问题求解未知量的具体表达式如下：
式中，
是第j个格网经过k+1次迭代后的电子密度变化值；P是第j个格网内经过的射线数；λ是松弛因子，用于控制收敛性，取值为0<λ<1；yi是第i条射线的dSTEC值；Ai,j是第i条射线在第j个格网中的截距；IEDj是第j个格网的电子密度IRI模型值；D是第i条射线上的修正量；W是第i条射线上第j个格网的dSTEC修正权重。
7.一种面向扰动探测的残差电离层层析方法装置，其特征在于，所述装置包括：
获取单元，被配置为获取观测站数据，所述观测站数据包括观测站的观测值文件、GPS导航电文文件、精密星历和差分码偏差文件，根据地面站的GPS观测数据计算STEC的观测值；
比较单元，被配置为基于穿刺点时间序列求电离层残差，采用区域电离层建模或GIM数据内插方法计算该穿刺点在同一时刻不同日期的VTEC，得到穿刺点VTEC时间序列后，通过多项式拟合计算研究时刻的拟合值，并通过投影函数得到STEC的拟合值，将STEC的拟合值与实测值比较得到STEC的异常变化；
建立单元，被配置为在层析区域建立格网模型，令每个体格网中的电子密度相同，建立残差电离层层析方程；
求解单元，被配置为利用同时代数重构技术解决层析方程组秩亏问题求解未知量。
8.根据权利要求7所述的装置，其特征在于，所述获取单元被进一步配置为通过GPS的双频信号计算得到所述STEC的观测值：
式中，f1和f2是GPS信号的频率，f1＝1575.42MHz,f2＝1227.60MHz,
和/>
是f1和f2频率信号的载波相位平滑伪距值；Δbk和Δbs分别是接收机和卫星的差分码偏差。
9.根据权利要求7所述的装置，其特征在于，所述比较单元被进一步配置为基于穿刺点时间序列求电离层残差，采用区域电离层建模或GIM数据内插方法计算该穿刺点在同一时刻不同日期的VTEC，得到穿刺点VTEC时间序列后，通过多项式拟合计算研究时刻的拟合值，并通过投影函数得到STEC的拟合值，将STEC的拟合值与实测值比较得到STEC的异常变化的具体表达式为：
式中，
为通过区域建模或GIM内插得到的ti时刻穿刺点VTEC值；VTEC(β,s)和VTEC(β,λ)分别是电离层区域模型和GIM内插，β为地理纬度，s为太阳时角，λ为地理经度，
为研究时刻t0的穿刺点VTEC拟合值；an为待求参数；/>
是穿刺点处STEC实测值；f(Z)是投影函数；/>
为求得穿刺点的STEC残差；Z为测站处卫星天顶角，Z′为穿刺点处卫星天顶角。
10.一种可读存储介质，其特征在于，所述可读存储介质存储有一个或者多个程序，所述一个或者多个程序可被一个或者多个处理器执行，以实现如权利要求1至7中任一项所述的方法。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310517850.8A CN116299595B (zh) | 2023-05-10 | 2023-05-10 | 面向扰动探测的残差电离层层析方法、装置及介质 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310517850.8A CN116299595B (zh) | 2023-05-10 | 2023-05-10 | 面向扰动探测的残差电离层层析方法、装置及介质 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN116299595A true CN116299595A (zh) | 2023-06-23 |
| CN116299595B CN116299595B (zh) | 2024-01-30 |
Family
ID=86788963
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310517850.8A Active CN116299595B (zh) | 2023-05-10 | 2023-05-10 | 面向扰动探测的残差电离层层析方法、装置及介质 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116299595B (zh) |
Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110275185A (zh) * | 2019-07-11 | 2019-09-24 | 武汉大学 | 基于gnss和geo卫星的电离层投影函数建模方法 |
| CN113093224A (zh) * | 2021-02-18 | 2021-07-09 | 北京航空航天大学 | 一种边缘增强的电离层层析方法 |
| CN113985454A (zh) * | 2021-10-23 | 2022-01-28 | 闽江学院 | 一种顾及方位角的电离层投影函数模型的建模方法 |
| CN114355419A (zh) * | 2021-12-15 | 2022-04-15 | 中国科学院国家授时中心 | 一种分布式北斗位置服务中心rtk产品定位方法及定位装置 |
| CN114355400A (zh) * | 2021-03-18 | 2022-04-15 | 武汉星辰北斗科技有限公司 | 一种基于电离层经向变化特征的电离层建模方法 |
| CN114690208A (zh) * | 2020-12-31 | 2022-07-01 | 复旦大学 | 电离层三维电子密度稀疏层析方法及其装置 |
| CN114721012A (zh) * | 2022-05-24 | 2022-07-08 | 武汉大学 | 一种电离层延迟效应改正方法及装置 |
-
2023
- 2023-05-10 CN CN202310517850.8A patent/CN116299595B/zh active Active
Patent Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110275185A (zh) * | 2019-07-11 | 2019-09-24 | 武汉大学 | 基于gnss和geo卫星的电离层投影函数建模方法 |
| CN114690208A (zh) * | 2020-12-31 | 2022-07-01 | 复旦大学 | 电离层三维电子密度稀疏层析方法及其装置 |
| CN113093224A (zh) * | 2021-02-18 | 2021-07-09 | 北京航空航天大学 | 一种边缘增强的电离层层析方法 |
| CN114355400A (zh) * | 2021-03-18 | 2022-04-15 | 武汉星辰北斗科技有限公司 | 一种基于电离层经向变化特征的电离层建模方法 |
| CN113985454A (zh) * | 2021-10-23 | 2022-01-28 | 闽江学院 | 一种顾及方位角的电离层投影函数模型的建模方法 |
| CN114355419A (zh) * | 2021-12-15 | 2022-04-15 | 中国科学院国家授时中心 | 一种分布式北斗位置服务中心rtk产品定位方法及定位装置 |
| CN114721012A (zh) * | 2022-05-24 | 2022-07-08 | 武汉大学 | 一种电离层延迟效应改正方法及装置 |
Non-Patent Citations (3)
| Title |
|---|
| CHANGZHI ZHAI 等: "Three‑dimensional reconstruction of seismo-traveling ionospheric disturbances after March 11, 2011, Japan Tohoku earthquake", JOURNAL OF GEODESY, pages 1 - 12 * |
| K. HEKI 等: "Mw dependence of the preseismic ionospheric electron enhancements", JOURNAL OF GEOPHYSICAL RESEARCH: SPACE PHYSICS, pages 7007 * |
| XULEI JIN 等: "Near real-time global ionospheric total electron content modeling and nowcasting based on GNSS observations", JOURNAL OF GEODESY, vol. 97, no. 3, pages 1 - 18 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN116299595B (zh) | 2024-01-30 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Rovira-Garcia et al. | A worldwide ionospheric model for fast precise point positioning | |
| EP3195011B1 (en) | Ionospheric scintillation prediction | |
| Yao et al. | An improved iterative algorithm for 3-D ionospheric tomography reconstruction | |
| CN110221320B (zh) | 一种基于抛物面天线观测的北斗频间偏差测定方法 | |
| Sharma et al. | Total solar eclipse of July 22, 2009: Its impact on the total electron content and ionospheric electron density in the Indian zone | |
| Yang et al. | Low-latitude GNSS ionospheric scintillation dependence on magnetic field orientation and impacts on positioning | |
| Komjathy et al. | An assessment of the current WAAS ionospheric correction algorithm in the south American region | |
| Seemala | Estimation of ionospheric total electron content (TEC) from GNSS observations | |
| Panda et al. | Ionospheric TEC variations at low latitude Indian region | |
| Hein et al. | Estimation method of ionospheric TEC distribution using single frequency measurements of GPS signals | |
| Huang et al. | Analysis and improvement of ionospheric thin shell model used in SBAS for China region | |
| Mohanty et al. | Ionospheric scintillation observation using space-borne synthetic aperture radar data | |
| Jiang et al. | Accuracy analysis of continuous deformation monitoring using BeiDou Navigation Satellite System at middle and high latitudes in China | |
| Radzi et al. | Seasonal variation of total electron content at equatorial station, Langkawi, Malaysia | |
| Jerez et al. | Two‐way assessment of ionospheric maps performance over the Brazilian region: Global versus regional products | |
| CN114384564A (zh) | 一种基于多源数据驱动的电离层层析成像方法 | |
| Yu et al. | Using the GPS observations to reconstruct the ionosphere three-dimensionally with an ionospheric data assimilation and analysis system (IDAAS) | |
| Bahadur et al. | Real-time single-frequency multi-GNSS positioning with ultra-rapid products | |
| Cheng et al. | Global monitoring of geomagnetic storm-induced ionosphere anomalies using 3-D ionospheric modeling with multi-GNSS and COSMIC measurements | |
| CN116299595B (zh) | 面向扰动探测的残差电离层层析方法、装置及介质 | |
| Cao et al. | Unveiling the driving factors of urban land subsidence in Beijing, China | |
| CN112083444B (zh) | 一种计及等离子体泡的低纬度机场电离层短时预测方法 | |
| Yu et al. | On the occurrence of F region irregularities over Haikou retrieved from COSMIC GPS radio occultation and ground-based ionospheric scintillation monitor observations | |
| Janssen | Likely impact of the approaching solar maximum on GNSS surveys: Be alert but not alarmed | |
| Geng et al. | The distribution characteristics of GPS cycle slip over the China mainland and adjacent region during the declining solar activity (2015–2018) period of solar cycle 24 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |