CN114384564A - 一种基于多源数据驱动的电离层层析成像方法 - Google Patents
一种基于多源数据驱动的电离层层析成像方法 Download PDFInfo
- Publication number
- CN114384564A CN114384564A CN202111678437.7A CN202111678437A CN114384564A CN 114384564 A CN114384564 A CN 114384564A CN 202111678437 A CN202111678437 A CN 202111678437A CN 114384564 A CN114384564 A CN 114384564A
- Authority
- CN
- China
- Prior art keywords
- ionosphere
- nmf2
- hmf2
- model
- data
- 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 112
- 238000003325 tomography Methods 0.000 title claims abstract description 46
- 238000000034 method Methods 0.000 title claims abstract description 22
- 238000012545 processing Methods 0.000 claims abstract description 8
- 238000007781 pre-processing Methods 0.000 claims abstract description 4
- 230000008859 change Effects 0.000 claims description 22
- 238000004364 calculation method Methods 0.000 claims description 18
- 230000000694 effects Effects 0.000 claims description 18
- 239000013256 coordination polymer Substances 0.000 claims description 15
- 239000011159 matrix material Substances 0.000 claims description 15
- 238000004422 calculation algorithm Methods 0.000 claims description 13
- 230000004907 flux Effects 0.000 claims description 12
- 238000005259 measurement Methods 0.000 claims description 9
- 238000009825 accumulation Methods 0.000 claims description 6
- JJWKPURADFRFRB-UHFFFAOYSA-N carbonyl sulfide Chemical compound O=C=S JJWKPURADFRFRB-UHFFFAOYSA-N 0.000 claims description 6
- 238000012821 model calculation Methods 0.000 claims description 6
- 238000009826 distribution Methods 0.000 claims description 4
- 238000000819 phase cycle Methods 0.000 claims description 4
- 230000005855 radiation Effects 0.000 claims description 4
- 230000002159 abnormal effect Effects 0.000 claims description 3
- 238000007405 data analysis Methods 0.000 claims description 3
- 230000001934 delay Effects 0.000 claims description 3
- 230000002452 interceptive effect Effects 0.000 claims description 3
- 230000005389 magnetism Effects 0.000 claims description 3
- 238000011160 research Methods 0.000 claims description 3
- 230000001186 cumulative effect Effects 0.000 claims description 2
- 230000003628 erosive effect Effects 0.000 claims description 2
- 230000005358 geomagnetic field Effects 0.000 claims description 2
- 238000009499 grossing Methods 0.000 claims description 2
- RAEKUJQUZFRYFZ-HVSMRNTNSA-N (4r,4ar,7r,7as)-3-methyl-7-[(5-nitropyridin-2-yl)disulfanyl]-2,4,4a,5,6,7,7a,13-octahydro-1h-4,12-methanobenzofuro[3,2-e]isoquinoline-9-ol Chemical compound S([C@@H]1CC[C@H]2[C@H]3CC=4C5=C(C(=CC=4)O)O[C@H]1C52CCN3C)SC1=CC=C([N+]([O-])=O)C=N1 RAEKUJQUZFRYFZ-HVSMRNTNSA-N 0.000 claims 3
- 238000013170 computed tomography imaging Methods 0.000 description 4
- 230000008569 process Effects 0.000 description 3
- 238000004891 communication Methods 0.000 description 2
- 125000004122 cyclic group Chemical group 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 238000012544 monitoring process Methods 0.000 description 2
- 238000010521 absorption reaction Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000005540 biological transmission Effects 0.000 description 1
- 238000004587 chromatography analysis Methods 0.000 description 1
- 238000012937 correction Methods 0.000 description 1
- 238000005260 corrosion Methods 0.000 description 1
- 230000007797 corrosion Effects 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 238000003384 imaging method Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 238000009828 non-uniform distribution Methods 0.000 description 1
- 230000001151 other effect Effects 0.000 description 1
- 230000035515 penetration Effects 0.000 description 1
- 230000000007 visual 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
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T11/00—2D [Two Dimensional] image generation
- G06T11/003—Reconstruction from projections, e.g. tomography
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T3/00—Geometric image transformations in the plane of the image
- G06T3/40—Scaling of whole images or parts thereof, e.g. expanding or contracting
- G06T3/4007—Scaling of whole images or parts thereof, e.g. expanding or contracting based on interpolation, e.g. bilinear interpolation
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Signal Processing (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Theoretical Computer Science (AREA)
- Computer Networks & Wireless Communication (AREA)
- Position Fixing By Use Of Radio Waves (AREA)
Abstract
本发明公开了一种基于多源数据驱动的电离层层析成像方法,包括如下步骤:步骤A,测高仪与掩星观测数据下载与预处理;步骤B,电离层F2层峰值电子密度NmF2驱动更新;步骤C,电离层F2层峰值高度hmF2驱动更新;步骤D,背景电离层模型电子密度驱动更新;步骤E,GNSS与卫星信标观测数据处理;步骤F,电离层电子密度层析成像。本发明所公开的方法,利用电离层测高仪和掩星垂直观测分辨率较高的特点,利用数据驱动的方法实现背景电离层模型的有效更新,在此基础上利用地基GNSS、卫星信标观测的TEC进行电离层层析成像,可以有效提升三维电离层电子密度的重构精度。
Description
技术领域
本发明属于电离层环境遥感技术领域,特别涉及该领域中的一种基于多源数据驱动的电离层层析成像方法。
背景技术
电离层是指地球上空60公里到1000公里高度的大气电离区域,穿越该区域的无线电波会产生折射、反射、散射和吸收等效应,从而影响导航、通信、测控、遥感等诸多无线电信息系统的性能。利用各类技术手段来获取电离层的特征参量,研究电离层中的各种现象,从而揭示其背后的物理机制,具有重要的科学乃至政治、军事和经济意义。电子密度是表征电离层状态最为直接的特征参量之一,获取精确的电子密度变化即可量化电离层对卫星导航、通信、雷达的影响。
电离层层析成像(电离层CT)是随着计算机硬件性能的提升和卫星电离层探测的进步而逐渐发展起来的一种非常重要的电离层遥感工具。电离层CT是利用卫星测量的电离层积分总电子含量(TEC)反演得到电子密度的过程。目前,可用于电离层层析的主要数据来源包括电离层测高仪(ionosonde)、地基LEO卫星信标(beacon)、地基GNSS、天基GNSS掩星(occultation)等数据。早期电离层CT信号来源主要是低轨信标卫星,随着GPS、GLONASS、GALILEO、北斗等全球卫星导航系统的发展,基于地基的GNSS台站逐步成为电离层CT主要数据来源。但受限于扫描视角和台站数目非均匀分布的影响,基于地基数据的电离层CT结果的质量极其依赖于所采用的背景电离层模型的精度,且不同程度存在着垂直分辨率较低、算法稳定性不足等问题。
发明内容
本发明所要解决的技术问题就是针对区域精细化电离层三维电子密度的监测需求,提供一种基于多源数据驱动的电离层层析成像方法。
本发明采用如下技术方案:
一种基于多源数据驱动的电离层层析成像方法,其改进之处在于,包括如下步骤:
步骤A,测高仪与掩星观测数据下载与预处理:
步骤A1,下载测高仪自动判读数据,数据来源包括全球无线电电离层观测站、空间物理交互式数据资源网和中国电波传播研究所测高仪观测网;
步骤A2,下载COSMIC数据分析和档案中心发布的COSMIC卫星、GRACE卫星、电子密度剖面数据产品;
步骤A3,提取全部观测数据的电离层F2峰值电子密度NmF2数值和峰值高度hmF2数值,分别标记为{x1,x2,...,xn}和{y1,y2,...,yn},其中n表示观测样本总数;
步骤A4,下载观测数据对应日期的太阳活动指数R12和电离层活动指数IG12;
步骤A5,将R12和IG12输入到国际参考电离层模型IRI中,计算得到对应时刻和对应位置处的F2峰值电子密度NmF2数据和峰值高度hmF2数据,分别标记为{xIRI,1,xIRI,2,...,xIRI,n}和{yIRI,1,yIRI,2,...,yIRI,n};
步骤A6,判定NmF2数据和峰值高度hmF2数据是否满足以下条件:
若满足以上条件,则保存该笔数据,否则认为数据存在较大偏差,剔除该笔数据;
步骤B,电离层F2层峰值电子密度NmF2驱动更新:
步骤B1,利用诺伊斯特雷利兹全球电离层峰值密度模型NPDM进行全球NmF2的变化特征建模,模型表述如下:
NmF2NPDM=F1·F2·F3·F4·F5
其中:F1项描述了电离层NmF2的逐日变化;F2项描述了电离层NmF2的年变化;F3描述了电离层NmF2随磁纬的变化;F4描述了电离层NmF2随磁赤道驼峰区的变化;F5描述了电离层NmF2随太阳活动的变化;
步骤B2,计算太阳赤经δ:
其中:ζ=180/π,DOY表示年积日;
步骤B3,计算电离层NmF2随地方时LT的时间变化项F1,F1描述了NmF2存在的日、半日和1/3日三个周期变化:
F1=P1+P2(c1cosVD+c2cosVSD+c3sinVSD+c4cosVTD+c5sinVTD)+c6BO
其中,LT表示地方时,VD、VSD、VTD分别表示三个不同谐波成分的角相位;P1和P2表示NmF2与太阳天顶角χ之间的依赖性;表示地理纬度,单位弧度;表示地理纬度,单位度;δ为太阳赤经,BO表示中纬区域电离层夏季咬蚀效应;LTBO表示咬蚀效应出现的地方时;c1—c6为模型系数;
步骤B4,计算电离层NmF2随年积日的时间变化项F2项,F2项描述了NmF2对年和半年的周期变化:
F2=1+c7cos(VA)+c8cos(VSA)
其中,VA表示年周期变化成分,VSA表示半年周期变化成分,c7—c8为模型系数;
步骤B5,计算电离层NmF2随磁纬的空间变化项F3项,F3描述了NmF2对地磁纬度的依赖性:
其中,和λ分别表示观测点的地理纬度和地理经度坐标,和分别表示北极磁极点的纬度和经度坐标,表示观测点的地磁纬度,c9为模型系数;地磁极点会随时间变化而变化,位置不是固定的常数,通过国际参考地磁场模型获得;
步骤B6,计算电离层NmF2随磁赤道驼峰区的空间变化项F4,利用高斯函数描述NmF2与磁赤道两侧的驼峰区位置之间的联系:
步骤B7,计算电离层NmF2随太阳活动的变化项F5,F5描述了NmF2与太阳射电通量F10.7指数之间的关联:
其中,F10.7表示太阳射电通量F10.7指数,c12—c13为模型系数;
步骤B8,利用垂测和掩星观测的NmF2,对诺伊斯特雷利兹全球电离层NmF2模型进行驱动更新,得到模型的拟合系数C(k)=(c1 (k),c2 (k),...,c13 (k)),计算方法如下:
fim(C(k))=|NmF2obs-NmF2NPDM(C(k))|im
其中,NmF2obs表示步骤A提取的测高仪和掩星NmF2观测值,NmF2NPDM表示NPDM模型计算的NmF2值,fim(C(k))表示数据驱动的极小化函数,下标im表示观测样本号,Num表示NmF2总的观测样本数,上标k表示迭代轮次,上标T表示矩阵转置;
步骤B9,将驱动更新后得到的拟合系数C重新代回步骤B1,利用NPDM模型计算生成层析成像对应区域全部网格点的NmF2值;
步骤C,电离层F2层峰值高度hmF2驱动更新:
步骤C1,利用诺伊斯特雷利兹全球电离层峰值高度模型NPHM进行全球hmF2的变化特征建模,模型表述如下:
hmF2NPHM=R1·R2·R3·R4
其中,R1项描述了电离层hmF2的逐日变化;R2项描述了电离层hmF2的年变化;R3描述了电离层hmF2随空间的变化;R4描述了电离层随太阳活动的变化;
步骤C2,计算太阳赤经δ:
其中:ζ=180/π,DOY表示年积日;
步骤C3,计算电离层hmF2随地方时LT的时间变化项R1,R1描述了hmF2存在的日、半日和1/3日三个周期变化:
R1=1+ρ1Z1+Z2(ρ2cosVD+ρ3cosVSD+ρ4sinVSD+ρ5cosVTD+ρ6sinVTD+ρ7cosVTD)
步骤C4:计算电离层hmF2随年积日的时间变化项R2项,R2项描述了hmF2对年和半年的周期变化:
R2=1+c8cos(VA)+c9cos(VSA)
其中,VA表示年周期变化成分,VSA表示半年周期变化成分,DOY表示年积日,ρ8—ρ9为模型系数;
步骤C5,计算电离层hmF2空间变化项R3,计算方法如下:
步骤C6,计算电离层hmF2随太阳活动的变化项R4,R4描述了hmF2与太阳射电通量F10.7指数之间的关联:
其中,F10.7表示太阳射电通量F10.7指数,ρ12—ρ13为模型系数;
步骤C7,利用垂测和掩星观测的hmF2,对NPHM模型进行驱动更新,得到模型的拟合系数CP(k)=(ρ1 (k),ρ2 (k),...,ρ13 (k)),计算方法如下:
fie(CP(k))=|hmF2obs-hmF2NPHM(CP(k))|ie
其中,hmF2obs表示步骤A提取的测高仪和掩星hmF2观测值,hmF2NPHM表示NPHM模型计算的hmF2值,fie(CP(k))表示数据驱动的极小化函数,下标ie表示观测样本号,Obs表示hmF2总的观测样本数,上标k表示迭代轮次,上标T表示矩阵转置;
步骤C8,将驱动更新后得到的拟合系数CP重新代回步骤C1,利用NPHM模型计算生成层析成像对应区域全部网格点的hmF2值;
步骤D,背景电离层模型电子密度驱动更新:
步骤D1,选择国际参考电离层模型IRI为电离层层析成像的背景电离层模型,设定模型输入的时间与步骤A1电离层测高仪和掩星接收机的观测时刻一致;
步骤D2,设定IRI模型的NmF2和hmF2为手动输入模式,分别将步骤B9和步骤C8计算得到的全部网格点的NmF2和hmF2值输入到IRI模型中;
步骤D3,利用驱动更新后的IRI模型计算并输出层析成像对应区域的三维电子密度值到文本文件中,输出到指定目录保存,此数据作为层析成像的背景电子密度值;
步骤E,GNSS与卫星信标观测数据处理:
步骤E1,对GNSS的观测数据进行周跳和异常值的探测,并采用双频的载波相位观测数据平滑码伪距,最后输出平滑处理后的观测数据;
步骤E2,采用Ciraolo L提出的方法,分别计算出卫星和接收机的硬件延迟,进而计算GNSS电离层倾斜TEC;
步骤E3,对卫星信标观测数据进行处理,读取差分多普勒相位序列,并获得卫星信标观测的电离层倾斜TEC;
步骤F,电离层电子密度层析成像:
步骤F1,构建电离层层析成像线性方程,方程具体表示如下:
式中,G是地基GNSS观测的射线数目,B是卫星信标观测的射线数目,N是总网格数,向量dd由地基GNSS和卫星信标接收机的TEC观测数据组成;AA是GNSS和卫星信标无线电信号传播路径在网格中的投影矩阵;X是未知电子密度的分布;e是测量与离散化后引入的误差;
步骤F2,对网格内同一高度层的电子密度加入水平约束,水平约束采用双线性插值算法实现:
其中,j表示网格编号,N为总网格数,w表示双线性插值所用的加权系数;
步骤F3,对步骤F1和F2中的线性方程组进行整合,同时忽略离散化误差,层析成像转换为求解下列线性方程组的问题:
d=AX
式中,向量d由dd和0值组合而成,A是投影矩阵AA和加权系数矩阵的组合;
步骤F4,读取步骤D3存储的电子密度数据文件,将其设定为电离层电子密度初值X(0);
步骤F5,利用代数重构算法ART进行电离层三维层析成像,ART算法计算方法如下:
其中,i为TEC路径编号;j为网格编号;和分别为迭代k+1次和k次的第j个网格的电子密度;aij为第i条路径在第j个网格内的投影;||ai||为第i条射线路径的总长;di为第i条射线路径的总电子含量;ξ为松弛因子;
步骤F6,ξk采用与迭代轮次k相关的函数,计算方法如下:
步骤F7,设置迭代截止门限tol,当tol<10-8时停止迭代,迭代截止门限计算方法为:
其中,N是总网格数;
步骤F8,ART算法迭代终止后,输出最终的X(k+1)到文件中并存储,即为最终的电离层电子密度层析成像结果。
本发明的有益效果是:
本发明所公开的方法,利用电离层测高仪和掩星垂直观测分辨率较高的特点,利用数据驱动的方法实现背景电离层模型的有效更新,在此基础上利用地基GNSS、卫星信标观测的TEC进行电离层层析成像,可以有效提升三维电离层电子密度的重构精度。本发明方法可为实现区域大范围高精度电离层监测提供技术支撑。
附图说明
图1是本发明方法的流程示意图。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白,以下结合附图和实施例,对本发明进行进一步详细说明。应当理解,此处所描述的具体实施例仅仅用以解释本发明,并不用于限定本发明。
在实际系统应用过程中,如卫星导航修正,短波无线电传播中,有时必须提供区域乃至全球三维分布的电离层电子密度信息,电离层三维CT技术就是实现这一目的的重要途径之一,电离层三维层析成像目前主要的数据来源是地基GNSS观测。由于GNSS卫星轨道高度较高的原因,对于某特定区域而言,绝大部分观测数据对应的射线路径相较于CT成像网格都是近垂直方向的,加上GNSS观测对于地球而言存在视角不完整的问题,按照电离层CT成像理论,这将导致地基GNSS成像所能达到的垂直分辨率较低,即反演得到的电离层峰值精度较差。为提高电离层三维CT成像的垂直分辨率,一个重要的解决方案是将多种地基和天基的观测数据源联合起来进行电离层CT成像。
实施例1,如图1所示,本实施例公开了一种基于多源数据驱动的电离层层析成像方法,包括如下步骤:
步骤A,测高仪与掩星观测数据下载与预处理:
步骤A1,下载测高仪自动判读数据,数据来源包括全球无线电电离层观测站(Global Ionospheric Radio Observatory,GIRO)、空间物理交互式数据资源网(TheSpace Physics Interactive Data Resource,SPIDR)和中国电波传播研究所(ChinaResearch Institute of Radiowave Propagation,CRIRP)测高仪观测网;
步骤A2,下载COSMIC数据分析和档案中心(COSMIC Data Analysis and ArchiveCenter,CDAAC)发布的COSMIC卫星、GRACE卫星、电子密度剖面数据产品(IonPrf);
步骤A3,提取全部观测数据的电离层F2峰值电子密度NmF2数值和峰值高度hmF2数值,分别标记为{x1,x2,...,xn}和{y1,y2,...,yn},其中n表示观测样本总数;
步骤A4,下载观测数据对应日期的太阳活动指数R12和电离层活动指数IG12;
步骤A5,将R12和IG12输入到国际参考电离层模型IRI中,计算得到对应时刻和对应位置处的F2峰值电子密度NmF2数据和峰值高度hmF2数据,分别标记为{xIRI,1,xIRI,2,...,xIRI,n}和{yIRI,1,yIRI,2,...,yIRI,n};
步骤A6,判定NmF2数据和峰值高度hmF2数据是否满足以下条件:
若满足以上条件,则保存该笔数据,否则认为数据存在较大偏差,剔除该笔数据;
步骤B,电离层F2层峰值电子密度NmF2驱动更新:
步骤B1,利用诺伊斯特雷利兹全球电离层峰值密度模型NPDM(Neustrelitz peakdensity model)进行全球NmF2的变化特征建模,模型表述如下:
NmF2NPDM=F1·F2·F3·F4·F5
其中:F1项描述了电离层NmF2的逐日变化;F2项描述了电离层NmF2的年变化;F3描述了电离层NmF2随磁纬的变化;F4描述了电离层NmF2随磁赤道驼峰区的变化;F5描述了电离层NmF2随太阳活动的变化;
步骤B2,计算太阳赤经δ:
其中:ζ=180/π,DOY表示年积日;
步骤B3,计算电离层NmF2随地方时LT的时间变化项F1,F1描述了NmF2存在的日、半日和1/3日三个周期变化:
F1=P1+P2(c1cosVD+c2cosVSD+c3sinVSD+c4cosVTD+c5sinVTD)+c6BO
其中,LT表示地方时,VD、VSD、VTD分别表示三个不同谐波成分的角相位;P1和P2表示NmF2与太阳天顶角χ之间的依赖性;表示地理纬度,单位弧度;表示地理纬度,单位度;δ为太阳赤经,BO表示中纬区域电离层夏季咬蚀(bite out)效应;LTBO表示咬蚀效应出现的地方时;c1—c6为模型系数;
步骤B4,计算电离层NmF2随年积日的时间变化项F2项,F2项描述了NmF2对年和半年的周期变化:
F2=1+c7cos(VA)+c8cos(VSA)
其中,VA表示年周期变化成分,VSA表示半年周期变化成分,c7—c8为模型系数;
步骤B5,计算电离层NmF2随磁纬的空间变化项F3项,F3描述了NmF2对地磁纬度的依赖性:
其中,和λ分别表示观测点的地理纬度和地理经度坐标,和分别表示北极磁极点的纬度和经度坐标,表示观测点的地磁纬度,c9为模型系数;地磁极点会随时间变化而变化,因此该位置不是固定的常数,可通过国际参考地磁场模型(InternationalGeomagnetic Reference Field model,IGRF)获得;
步骤B6,计算电离层NmF2随磁赤道驼峰区的空间变化项F4,利用高斯函数描述NmF2与磁赤道两侧的驼峰区位置之间的联系:
步骤B7,计算电离层NmF2随太阳活动的变化项F5,F5描述了NmF2与太阳射电通量F10.7指数之间的关联:
其中,F10.7表示太阳射电通量F10.7指数,c12—c13为模型系数;
步骤B8,利用垂测和掩星观测的NmF2,对诺伊斯特雷利兹全球电离层NmF2模型进行驱动更新,得到模型的拟合系数C(k)=(c1 (k),c2 (k),...,c13 (k)),计算方法如下:
fim(C(k))=|NmF2obs-NmF2NPDM(C(k))|im
其中,NmF2obs表示步骤A提取的测高仪和掩星NmF2观测值,NmF2NPDM表示NPDM模型计算的NmF2值,fim(C(k))表示数据驱动的极小化函数,下标im表示观测样本号,Num表示NmF2总的观测样本数,上标k表示迭代轮次,上标T表示矩阵转置;
步骤B9,将驱动更新后得到的拟合系数C重新代回步骤B1,利用NPDM模型计算生成层析成像对应区域全部网格点的NmF2值;
步骤C,电离层F2层峰值高度hmF2驱动更新:
步骤C1,利用诺伊斯特雷利兹全球电离层峰值高度模型NPHM(Neustrelitz peakheight model)进行全球hmF2的变化特征建模,模型表述如下:
hmF2NPHM=R1·R2·R3·R4
其中,R1项描述了电离层hmF2的逐日变化;R2项描述了电离层hmF2的年变化;R3描述了电离层hmF2随空间的变化;R4描述了电离层随太阳活动的变化;
步骤C2,计算太阳赤经δ:
其中:ζ=180/π,DOY表示年积日;
步骤C3,计算电离层hmF2随地方时LT的时间变化项R1,R1描述了hmF2存在的日、半日和1/3日三个周期变化:
R1=1+ρ1Z1+Z2(ρ2cosVD+ρ3cosVSD+ρ4sinVSD+ρ5cosVTD+ρ6sinVTD+ρ7cosVTD)
步骤C4:计算电离层hmF2随年积日的时间变化项R2项,R2项描述了hmF2对年和半年的周期变化:
R2=1+c8cos(VA)+c9cos(VSA)
其中,VA表示年周期变化成分,VSA表示半年周期变化成分,DOY表示年积日,ρ8—ρ9为模型系数;
步骤C5,计算电离层hmF2空间变化项R3,计算方法如下:
步骤C6,计算电离层hmF2随太阳活动的变化项R4,R4描述了hmF2与太阳射电通量F10.7指数之间的关联:
其中,F10.7表示太阳射电通量F10.7指数,ρ12—ρ13为模型系数;
步骤C7,利用垂测和掩星观测的hmF2,对NPHM模型进行驱动更新,得到模型的拟合系数CP(k)=(ρ1 (k),ρ2 (k),...,ρ13 (k)),计算方法如下:
fie(CP(k))=|hmF2obs-hmF2NPHM(CP(k))|ie
其中,hmF2obs表示步骤A提取的测高仪和掩星hmF2观测值,hmF2NPHM表示NPHM模型计算的hmF2值,fie(CP(k))表示数据驱动的极小化函数,下标ie表示观测样本号,Obs表示hmF2总的观测样本数,上标k表示迭代轮次,上标T表示矩阵转置;
步骤C8,将驱动更新后得到的拟合系数CP重新代回步骤C1,利用NPHM模型计算生成层析成像对应区域全部网格点的hmF2值;
步骤D,背景电离层模型电子密度驱动更新:
步骤D1,选择国际参考电离层模型IRI为电离层层析成像的背景电离层模型,设定模型输入的时间与步骤A1电离层测高仪和掩星接收机的观测时刻一致;
步骤D2,设定IRI模型的NmF2和hmF2为手动输入模式,分别将步骤B9和步骤C8计算得到的全部网格点的NmF2和hmF2值输入到IRI模型中;
步骤D3,利用驱动更新后的IRI模型计算并输出层析成像对应区域的三维电子密度值到文本文件中,输出到指定目录保存,此数据作为层析成像的背景电子密度值;
步骤E,GNSS与卫星信标观测数据处理:
步骤E1,对GNSS观测数据进行处理,采用Bernese软件对GNSS的观测数据进行周跳和异常值的探测,并采用双频的载波相位观测数据平滑码伪距,最后平滑处理后的观测数据仍以RINEX格式的文件输出;
步骤E2,采用Ciraolo L提出的方法(具体参考Ciraolo L,Azpilicueta F,Brunini C,et al.Calibration errors on experimental slant total electroncontent(TEC)determined with GPS[J].J.Geod.,2007,81(2):111–120),分别计算出卫星和接收机的硬件延迟,进而计算GNSS电离层倾斜TEC;
步骤E3,提取三频信标接收机VHF、UHF和L频段同相支路I路和正交支路Q路数据、卫星两行轨道根数TLE星历数据;根据TLE星历计算卫星位置,同时记录接收机数据采集时刻;利用提取的接收机I、Q路数据求差分多普勒相位,并将差分多普勒相位与观测时刻、卫星位置进行匹配;读取差分多普勒相位序列,利用相位连接技术修复相位翻转,获取相位连接后的差分多普勒相位序列,计算得到电离层相对TEC;利用多站法计算未知相位积分常数,计算得到电离层倾斜绝对TEC;
步骤F,电离层电子密度层析成像:
步骤F1,构建电离层层析成像线性方程,方程具体表示如下:
式中,G是地基GNSS观测的射线数目,B是卫星信标观测的射线数目,N是总网格数,向量dd由地基GNSS和卫星信标接收机的TEC观测数据组成;AA是GNSS和卫星信标无线电信号传播路径在网格中的投影矩阵;X是未知电子密度的分布;e是测量与离散化后引入的误差;
步骤F2,对网格内同一高度层的电子密度加入水平约束,以实现对无射线穿越网格内电子密度的调整。水平约束采用双线性插值算法实现:
其中,j表示网格编号,N为总网格数,w表示双线性插值所用的加权系数;
步骤F3,对步骤F1和F2中的线性方程组进行整合,同时忽略离散化误差,层析成像转换为求解下列线性方程组的问题:
d=AX
式中,向量d由dd和0值组合而成,A是投影矩阵AA和加权系数矩阵的组合;
步骤F4,读取步骤D3存储的电子密度数据文件,将其设定为电离层电子密度初值X(0);
步骤F5,利用代数重构算法ART进行电离层三维层析成像,ART算法计算方法如下:
其中,i为TEC路径编号;j为网格编号;和分别为迭代k+1次和k次的第j个网格的电子密度;aij为第i条路径在第j个网格内的投影;||ai||为第i条射线路径的总长;di为第i条射线路径的总电子含量;ξ为松弛因子;
步骤F6,ξk采用与迭代轮次k相关的函数,计算方法如下:
步骤F7,设置迭代截止门限tol,当tol<10-8时停止迭代,迭代截止门限计算方法为:
其中,N是总网格数;
步骤F8,ART算法迭代终止后,输出最终的X(k+1)到文件中并存储,即为最终的电离层电子密度层析成像结果。
Claims (1)
1.一种基于多源数据驱动的电离层层析成像方法,其特征在于,包括如下步骤:
步骤A,测高仪与掩星观测数据下载与预处理:
步骤A1,下载测高仪自动判读数据,数据来源包括全球无线电电离层观测站、空间物理交互式数据资源网和中国电波传播研究所测高仪观测网;
步骤A2,下载COSMIC数据分析和档案中心发布的COSMIC卫星、GRACE卫星、电子密度剖面数据产品;
步骤A3,提取全部观测数据的电离层F2峰值电子密度NmF2数值和峰值高度hmF2数值,分别标记为{x1,x2,...,xn}和{y1,y2,...,yn},其中n表示观测样本总数;
步骤A4,下载观测数据对应日期的太阳活动指数R12和电离层活动指数IG12;
步骤A5,将R12和IG12输入到国际参考电离层模型IRI中,计算得到对应时刻和对应位置处的F2峰值电子密度NmF2数据和峰值高度hmF2数据,分别标记为{xIRI,1,xIRI,2,...,xIRI,n}和{yIRI,1,yIRI,2,...,yIRI,n};
步骤A6,判定NmF2数据和峰值高度hmF2数据是否满足以下条件:
若满足以上条件,则保存该笔数据,否则认为数据存在较大偏差,剔除该笔数据;
步骤B,电离层F2层峰值电子密度NmF2驱动更新:
步骤B1,利用诺伊斯特雷利兹全球电离层峰值密度模型NPDM进行全球NmF2的变化特征建模,模型表述如下:
NmF2NPDM=F1·F2·F3·F4·F5
其中:F1项描述了电离层NmF2的逐日变化;F2项描述了电离层NmF2的年变化;F3描述了电离层NmF2随磁纬的变化;F4描述了电离层NmF2随磁赤道驼峰区的变化;F5描述了电离层NmF2随太阳活动的变化;
步骤B2,计算太阳赤经δ:
其中:ζ=180/π,DOY表示年积日;
步骤B3,计算电离层NmF2随地方时LT的时间变化项F1,F1描述了NmF2存在的日、半日和1/3日三个周期变化:
F1=P1+P2(c1cosVD+c2cosVSD+c3sinVSD+c4cosVTD+c5sinVTD)+c6BO
其中,LT表示地方时,VD、VSD、VTD分别表示三个不同谐波成分的角相位;P1和P2表示NmF2与太阳天顶角χ之间的依赖性;表示地理纬度,单位弧度;表示地理纬度,单位度;δ为太阳赤经,BO表示中纬区域电离层夏季咬蚀效应;LTBO表示咬蚀效应出现的地方时;c1—c6为模型系数;
步骤B4,计算电离层NmF2随年积日的时间变化项F2项,F2项描述了NmF2对年和半年的周期变化:
F2=1+c7cos(VA)+c8cos(VSA)
其中,VA表示年周期变化成分,VSA表示半年周期变化成分,c7—c8为模型系数;
步骤B5,计算电离层NmF2随磁纬的空间变化项F3项,F3描述了NmF2对地磁纬度的依赖性:
其中,和λ分别表示观测点的地理纬度和地理经度坐标,和分别表示北极磁极点的纬度和经度坐标,表示观测点的地磁纬度,c9为模型系数;地磁极点会随时间变化而变化,位置不是固定的常数,通过国际参考地磁场模型获得;
步骤B6,计算电离层NmF2随磁赤道驼峰区的空间变化项F4,利用高斯函数描述NmF2与磁赤道两侧的驼峰区位置之间的联系:
步骤B7,计算电离层NmF2随太阳活动的变化项F5,F5描述了NmF2与太阳射电通量F10.7指数之间的关联:
其中,F10.7表示太阳射电通量F10.7指数,c12—c13为模型系数;
步骤B8,利用垂测和掩星观测的NmF2,对诺伊斯特雷利兹全球电离层NmF2模型进行驱动更新,得到模型的拟合系数C(k)=(c1 (k),c2 (k),...,c13 (k)),计算方法如下:
fim(C(k))=|NmF2obs-NmF2NPDM(C(k))|im
其中,NmF2obs表示步骤A提取的测高仪和掩星NmF2观测值,NmF2NPDM表示NPDM模型计算的NmF2值,fim(C(k))表示数据驱动的极小化函数,下标im表示观测样本号,Num表示NmF2总的观测样本数,上标k表示迭代轮次,上标T表示矩阵转置;
步骤B9,将驱动更新后得到的拟合系数C重新代回步骤B1,利用NPDM模型计算生成层析成像对应区域全部网格点的NmF2值;
步骤C,电离层F2层峰值高度hmF2驱动更新:
步骤C1,利用诺伊斯特雷利兹全球电离层峰值高度模型NPHM进行全球hmF2的变化特征建模,模型表述如下:
hmF2NPHM=R1·R2·R3·R4
其中,R1项描述了电离层hmF2的逐日变化;R2项描述了电离层hmF2的年变化;R3描述了电离层hmF2随空间的变化;R4描述了电离层随太阳活动的变化;
步骤C2,计算太阳赤经δ:
其中:ζ=180/π,DOY表示年积日;
步骤C3,计算电离层hmF2随地方时LT的时间变化项R1,R1描述了hmF2存在的日、半日和1/3日三个周期变化:
R1=1+ρ1Z1+Z2(ρ2cosVD+ρ3cosVSD+ρ4sinVSD+ρ5cosVTD+ρ6sinVTD+ρ7cosVTD)
步骤C4:计算电离层hmF2随年积日的时间变化项R2项,R2项描述了hmF2对年和半年的周期变化:
R2=1+c8cos(VA)+c9cos(VSA)
其中,VA表示年周期变化成分,VSA表示半年周期变化成分,DOY表示年积日,ρ8—ρ9为模型系数;
步骤C5,计算电离层hmF2空间变化项R3,计算方法如下:
步骤C6,计算电离层hmF2随太阳活动的变化项R4,R4描述了hmF2与太阳射电通量F10.7指数之间的关联:
其中,F10.7表示太阳射电通量F10.7指数,ρ12—ρ13为模型系数;
步骤C7,利用垂测和掩星观测的hmF2,对NPHM模型进行驱动更新,得到模型的拟合系数CP(k)=(ρ1 (k),ρ2 (k),...,ρ13 (k)),计算方法如下:
fie(CP(k))=|hmF2obs-hmF2NPHM(CP(k))|ie
其中,hmF2obs表示步骤A提取的测高仪和掩星hmF2观测值,hmF2NPHM表示NPHM模型计算的hmF2值,fie(CP(k))表示数据驱动的极小化函数,下标ie表示观测样本号,Obs表示hmF2总的观测样本数,上标k表示迭代轮次,上标T表示矩阵转置;
步骤C8,将驱动更新后得到的拟合系数CP重新代回步骤C1,利用NPHM模型计算生成层析成像对应区域全部网格点的hmF2值;
步骤D,背景电离层模型电子密度驱动更新:
步骤D1,选择国际参考电离层模型IRI为电离层层析成像的背景电离层模型,设定模型输入的时间与步骤A1电离层测高仪和掩星接收机的观测时刻一致;
步骤D2,设定IRI模型的NmF2和hmF2为手动输入模式,分别将步骤B9和步骤C8计算得到的全部网格点的NmF2和hmF2值输入到IRI模型中;
步骤D3,利用驱动更新后的IRI模型计算并输出层析成像对应区域的三维电子密度值到文本文件中,输出到指定目录保存,此数据作为层析成像的背景电子密度值;
步骤E,GNSS与卫星信标观测数据处理:
步骤E1,对GNSS的观测数据进行周跳和异常值的探测,并采用双频的载波相位观测数据平滑码伪距,最后输出平滑处理后的观测数据;
步骤E2,采用Ciraolo L提出的方法,分别计算出卫星和接收机的硬件延迟,进而计算GNSS电离层倾斜TEC;
步骤E3,对卫星信标观测数据进行处理,读取差分多普勒相位序列,并获得卫星信标观测的电离层倾斜TEC;
步骤F,电离层电子密度层析成像:
步骤F1,构建电离层层析成像线性方程,方程具体表示如下:
式中,G是地基GNSS观测的射线数目,B是卫星信标观测的射线数目,N是总网格数,向量dd由地基GNSS和卫星信标接收机的TEC观测数据组成;AA是GNSS和卫星信标无线电信号传播路径在网格中的投影矩阵;X是待求解的电子密度分布;e是测量与离散化后引入的误差;
步骤F2,对网格内同一高度层的电子密度加入水平约束,水平约束采用双线性插值算法实现:
其中,j表示网格编号,N为总网格数,w表示双线性插值所用的加权系数;
步骤F3,对步骤F1和F2中的线性方程组进行整合,同时忽略离散化误差,层析成像转换为求解下列线性方程组的问题:
d=AX
式中,向量d由dd和0值组合而成,A是投影矩阵AA和加权系数矩阵的组合;
步骤F4,读取步骤D3存储的电子密度数据文件,将其设定为电离层电子密度初值X(0);
步骤F5,利用代数重构算法ART进行电离层三维层析成像,ART算法计算方法如下:
其中,i为TEC路径编号;j为网格编号;和分别为迭代k+1次和k次的第j个网格的电子密度;aij为第i条路径在第j个网格内的投影;||ai||为第i条射线路径的总长;di为第i条射线路径的总电子含量;ξ为松弛因子;
步骤F6,ξk采用与迭代轮次k相关的函数,计算方法如下:
步骤F7,设置迭代截止门限tol,当tol<10-8时停止迭代,迭代截止门限计算方法为:
其中,N是总网格数;
步骤F8,ART算法迭代终止后,输出最终的X(k+1)到文件中并存储,即为最终的电离层电子密度层析成像结果。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202111678437.7A CN114384564B (zh) | 2021-12-31 | 2021-12-31 | 一种基于多源数据驱动的电离层层析成像方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202111678437.7A CN114384564B (zh) | 2021-12-31 | 2021-12-31 | 一种基于多源数据驱动的电离层层析成像方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN114384564A true CN114384564A (zh) | 2022-04-22 |
CN114384564B CN114384564B (zh) | 2024-05-14 |
Family
ID=81200116
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202111678437.7A Active CN114384564B (zh) | 2021-12-31 | 2021-12-31 | 一种基于多源数据驱动的电离层层析成像方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN114384564B (zh) |
Cited By (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN115078848A (zh) * | 2022-07-12 | 2022-09-20 | 武汉大学 | 基于闪电信号的电离层无源被动探测方法 |
CN115792963A (zh) * | 2023-02-13 | 2023-03-14 | 天津云遥宇航科技有限公司 | 基于gim模型的电离层校正方法、电子设备及存储介质 |
Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN111125609A (zh) * | 2019-12-20 | 2020-05-08 | 中国电波传播研究所(中国电子科技集团公司第二十二研究所) | 一种基于双指数驱动的电离层三维电子密度重构方法 |
CN112526617A (zh) * | 2020-11-19 | 2021-03-19 | 中国电波传播研究所(中国电子科技集团公司第二十二研究所) | 一种基于多源卫星信号的电离层层析成像系统观测数据模拟方法 |
US20210389472A1 (en) * | 2020-06-15 | 2021-12-16 | East China Jiaotong University | Computerized ionospheric tomography method based on vertical boundary truncation rays |
-
2021
- 2021-12-31 CN CN202111678437.7A patent/CN114384564B/zh active Active
Patent Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN111125609A (zh) * | 2019-12-20 | 2020-05-08 | 中国电波传播研究所(中国电子科技集团公司第二十二研究所) | 一种基于双指数驱动的电离层三维电子密度重构方法 |
US20210389472A1 (en) * | 2020-06-15 | 2021-12-16 | East China Jiaotong University | Computerized ionospheric tomography method based on vertical boundary truncation rays |
CN112526617A (zh) * | 2020-11-19 | 2021-03-19 | 中国电波传播研究所(中国电子科技集团公司第二十二研究所) | 一种基于多源卫星信号的电离层层析成像系统观测数据模拟方法 |
Non-Patent Citations (1)
Title |
---|
欧明;甄卫民;徐继生;刘钝;於晓;: "地基GPS与掩星联合的电离层层析成像方法研究", 全球定位系统, no. 05, 31 October 2014 (2014-10-31) * |
Cited By (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN115078848A (zh) * | 2022-07-12 | 2022-09-20 | 武汉大学 | 基于闪电信号的电离层无源被动探测方法 |
CN115078848B (zh) * | 2022-07-12 | 2024-04-09 | 武汉大学 | 基于闪电信号的电离层无源被动探测方法 |
CN115792963A (zh) * | 2023-02-13 | 2023-03-14 | 天津云遥宇航科技有限公司 | 基于gim模型的电离层校正方法、电子设备及存储介质 |
Also Published As
Publication number | Publication date |
---|---|
CN114384564B (zh) | 2024-05-14 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN111273335B (zh) | 一种基于垂测数据约束的电离层层析成像方法 | |
Yue et al. | Characterizing GPS radio occultation loss of lock due to ionospheric weather | |
Hajj et al. | COSMIC GPS ionospheric sensing and space weather | |
Hadas et al. | Impact and implementation of higher‐order ionospheric effects on precise GNSS applications | |
CN101403790A (zh) | 单频gps接收机的精密单点定位方法 | |
Pireaux et al. | Higher-order ionospheric effects in GPS time and frequency transfer | |
CN114384564B (zh) | 一种基于多源数据驱动的电离层层析成像方法 | |
Weiss et al. | Orbit and clock product generation | |
Seemala | Estimation of ionospheric total electron content (TEC) from GNSS observations | |
Prol et al. | Global‐scale ionospheric tomography during the March 17, 2015 geomagnetic storm | |
Ren et al. | Electron density reconstruction by ionospheric tomography from the combination of GNSS and upcoming LEO constellations | |
Martin | GNSS precise point positioning: The enhancement with GLONASS | |
Jakowski | Ionosphere monitoring | |
Conrad et al. | Improved GPS-based single-frequency orbit determination for the CYGNSS spacecraft using GipsyX | |
Muradyan et al. | GPS/INS navigation precision and its effect on airborne radio occultation retrieval accuracy | |
Jin et al. | Near real-time global ionospheric total electron content modeling and nowcasting based on GNSS observations | |
Syndergaard | A new algorithm for retrieving GPS radio occultation total electron content | |
Glaner | Towards instantaneous PPP convergence using multiple GNSS signals | |
Krynski et al. | National report of Poland to EUREF 2013 | |
Bourne | An algorithm for accurate ionospheric total electron content and receiver bias estimation using GPS measurements | |
Zhao et al. | Error analysis for the baseline estimation and calibration of distributed InSAR satellites | |
Nikoukar et al. | A novel data assimilation technique for the plasmasphere | |
Christovam et al. | PPP at low latitudes with ionospheric model exclusively based on single frequency GNSS measurements | |
Nohutcu | Development of a MATLAB based software package for ionosphere modeling | |
Wu et al. | Evaluation of Abel inversion method assisted by an improved IRI model |
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 |