CN114384564A - 一种基于多源数据驱动的电离层层析成像方法 - Google Patents

一种基于多源数据驱动的电离层层析成像方法 Download PDF

Info

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
Application number
CN202111678437.7A
Other languages
English (en)
Other versions
CN114384564B (zh
Inventor
欧明
甄卫民
冯健
吴家燕
陈龙江
王妍
许娜
陈亮
马宝田
刘钝
於晓
陈丽
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China Institute of Radio Wave Propagation CETC 22 Research Institute
Original Assignee
China Institute of Radio Wave Propagation CETC 22 Research Institute
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by China Institute of Radio Wave Propagation CETC 22 Research Institute filed Critical China Institute of Radio Wave Propagation CETC 22 Research Institute
Priority to CN202111678437.7A priority Critical patent/CN114384564B/zh
Publication of CN114384564A publication Critical patent/CN114384564A/zh
Application granted granted Critical
Publication of CN114384564B publication Critical patent/CN114384564B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01SRADIO 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/00Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/01Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/13Receivers
    • G01S19/35Constructional details or hardware or software details of the signal processing chain
    • G01S19/37Hardware or software details of the signal processing chain
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T11/002D [Two Dimensional] image generation
    • G06T11/003Reconstruction from projections, e.g. tomography
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T3/00Geometric image transformations in the plane of the image
    • G06T3/40Scaling of whole images or parts thereof, e.g. expanding or contracting
    • G06T3/4007Scaling 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数据是否满足以下条件:
Figure BDA0003453142400000021
Figure BDA0003453142400000022
若满足以上条件,则保存该笔数据,否则认为数据存在较大偏差,剔除该笔数据;
步骤B,电离层F2层峰值电子密度NmF2驱动更新:
步骤B1,利用诺伊斯特雷利兹全球电离层峰值密度模型NPDM进行全球NmF2的变化特征建模,模型表述如下:
NmF2NPDM=F1·F2·F3·F4·F5
其中:F1项描述了电离层NmF2的逐日变化;F2项描述了电离层NmF2的年变化;F3描述了电离层NmF2随磁纬的变化;F4描述了电离层NmF2随磁赤道驼峰区的变化;F5描述了电离层NmF2随太阳活动的变化;
步骤B2,计算太阳赤经δ:
Figure BDA0003453142400000023
其中:ζ=180/π,DOY表示年积日;
步骤B3,计算电离层NmF2随地方时LT的时间变化项F1,F1描述了NmF2存在的日、半日和1/3日三个周期变化:
F1=P1+P2(c1cosVD+c2cosVSD+c3sinVSD+c4cosVTD+c5sinVTD)+c6BO
Figure BDA0003453142400000031
Figure BDA0003453142400000032
Figure BDA0003453142400000033
Figure BDA0003453142400000034
Figure BDA0003453142400000035
其中,LT表示地方时,VD、VSD、VTD分别表示三个不同谐波成分的角相位;P1和P2表示NmF2与太阳天顶角χ之间的依赖性;
Figure BDA0003453142400000036
表示地理纬度,单位弧度;
Figure BDA0003453142400000037
表示地理纬度,单位度;δ为太阳赤经,BO表示中纬区域电离层夏季咬蚀效应;LTBO表示咬蚀效应出现的地方时;c1—c6为模型系数;
步骤B4,计算电离层NmF2随年积日的时间变化项F2项,F2项描述了NmF2对年和半年的周期变化:
F2=1+c7cos(VA)+c8cos(VSA)
Figure BDA0003453142400000038
Figure BDA0003453142400000039
其中,VA表示年周期变化成分,VSA表示半年周期变化成分,c7—c8为模型系数;
步骤B5,计算电离层NmF2随磁纬的空间变化项F3项,F3描述了NmF2对地磁纬度的依赖性:
Figure BDA00034531424000000310
Figure BDA00034531424000000311
其中,
Figure BDA0003453142400000041
和λ分别表示观测点的地理纬度和地理经度坐标,
Figure BDA0003453142400000042
Figure BDA0003453142400000043
分别表示北极磁极点的纬度和经度坐标,
Figure BDA0003453142400000044
表示观测点的地磁纬度,c9为模型系数;地磁极点会随时间变化而变化,位置不是固定的常数,通过国际参考地磁场模型获得;
步骤B6,计算电离层NmF2随磁赤道驼峰区的空间变化项F4,利用高斯函数描述NmF2与磁赤道两侧的驼峰区位置之间的联系:
Figure BDA0003453142400000045
Figure BDA0003453142400000046
其中:
Figure BDA0003453142400000047
Figure BDA0003453142400000048
分别表示北侧和南侧驼峰区磁纬度,分别设置为16°N和15°S,σc表示高斯半峰宽,LT表示地方时,c10—c11为模型系数;
步骤B7,计算电离层NmF2随太阳活动的变化项F5,F5描述了NmF2与太阳射电通量F10.7指数之间的关联:
Figure BDA0003453142400000049
其中,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
Figure BDA00034531424000000410
Figure BDA0003453142400000051
其中,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,计算太阳赤经δ:
Figure BDA0003453142400000052
其中:ζ=180/π,DOY表示年积日;
步骤C3,计算电离层hmF2随地方时LT的时间变化项R1,R1描述了hmF2存在的日、半日和1/3日三个周期变化:
R1=1+ρ1Z1+Z22cosVD3cosVSD4sinVSD5cosVTD6sinVTD7cosVTD)
Figure BDA0003453142400000061
Figure BDA0003453142400000062
Figure BDA0003453142400000063
其中,LT表示地方时,VD、VSD、VTD分别表示三个谐波成分的角相位;Z1和Z2表示hmF2与太阳天顶角χ之间的依赖性,其中
Figure BDA0003453142400000064
表示地理纬度,δ为太阳赤经,ρ1—ρ7为模型系数;
步骤C4:计算电离层hmF2随年积日的时间变化项R2项,R2项描述了hmF2对年和半年的周期变化:
R2=1+c8cos(VA)+c9cos(VSA)
Figure BDA0003453142400000065
Figure BDA0003453142400000066
其中,VA表示年周期变化成分,VSA表示半年周期变化成分,DOY表示年积日,ρ8—ρ9为模型系数;
步骤C5,计算电离层hmF2空间变化项R3,计算方法如下:
Figure BDA0003453142400000067
Figure BDA0003453142400000068
其中,σc表示高斯半峰宽,
Figure BDA0003453142400000069
Figure BDA00034531424000000610
分别设置为40°和20°,LT表示地方时,ρ10—ρ11为模型系数,σLT表示当地时14点对应的高斯半峰宽;
步骤C6,计算电离层hmF2随太阳活动的变化项R4,R4描述了hmF2与太阳射电通量F10.7指数之间的关联:
Figure BDA0003453142400000071
其中,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
Figure BDA0003453142400000072
Figure BDA0003453142400000073
其中,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,构建电离层层析成像线性方程,方程具体表示如下:
Figure BDA0003453142400000081
式中,G是地基GNSS观测的射线数目,B是卫星信标观测的射线数目,N是总网格数,向量dd由地基GNSS和卫星信标接收机的TEC观测数据组成;AA是GNSS和卫星信标无线电信号传播路径在网格中的投影矩阵;X是未知电子密度的分布;e是测量与离散化后引入的误差;
步骤F2,对网格内同一高度层的电子密度加入水平约束,水平约束采用双线性插值算法实现:
Figure BDA0003453142400000082
其中,j表示网格编号,N为总网格数,w表示双线性插值所用的加权系数;
步骤F3,对步骤F1和F2中的线性方程组进行整合,同时忽略离散化误差,层析成像转换为求解下列线性方程组的问题:
d=AX
式中,向量d由dd和0值组合而成,A是投影矩阵AA和加权系数矩阵的组合;
步骤F4,读取步骤D3存储的电子密度数据文件,将其设定为电离层电子密度初值X(0)
步骤F5,利用代数重构算法ART进行电离层三维层析成像,ART算法计算方法如下:
Figure BDA0003453142400000091
其中,i为TEC路径编号;j为网格编号;
Figure BDA0003453142400000092
Figure BDA0003453142400000093
分别为迭代k+1次和k次的第j个网格的电子密度;aij为第i条路径在第j个网格内的投影;||ai||为第i条射线路径的总长;di为第i条射线路径的总电子含量;ξ为松弛因子;
步骤F6,ξk采用与迭代轮次k相关的函数,计算方法如下:
Figure BDA0003453142400000094
步骤F7,设置迭代截止门限tol,当tol<10-8时停止迭代,迭代截止门限计算方法为:
Figure BDA0003453142400000095
其中,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数据是否满足以下条件:
Figure BDA0003453142400000101
Figure BDA0003453142400000102
若满足以上条件,则保存该笔数据,否则认为数据存在较大偏差,剔除该笔数据;
步骤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,计算太阳赤经δ:
Figure BDA0003453142400000111
其中:ζ=180/π,DOY表示年积日;
步骤B3,计算电离层NmF2随地方时LT的时间变化项F1,F1描述了NmF2存在的日、半日和1/3日三个周期变化:
F1=P1+P2(c1cosVD+c2cosVSD+c3sinVSD+c4cosVTD+c5sinVTD)+c6BO
Figure BDA0003453142400000112
Figure BDA0003453142400000113
Figure BDA0003453142400000114
Figure BDA0003453142400000115
Figure BDA0003453142400000116
其中,LT表示地方时,VD、VSD、VTD分别表示三个不同谐波成分的角相位;P1和P2表示NmF2与太阳天顶角χ之间的依赖性;
Figure BDA0003453142400000121
表示地理纬度,单位弧度;
Figure BDA0003453142400000122
表示地理纬度,单位度;δ为太阳赤经,BO表示中纬区域电离层夏季咬蚀(bite out)效应;LTBO表示咬蚀效应出现的地方时;c1—c6为模型系数;
步骤B4,计算电离层NmF2随年积日的时间变化项F2项,F2项描述了NmF2对年和半年的周期变化:
F2=1+c7cos(VA)+c8cos(VSA)
Figure BDA0003453142400000123
Figure BDA0003453142400000124
其中,VA表示年周期变化成分,VSA表示半年周期变化成分,c7—c8为模型系数;
步骤B5,计算电离层NmF2随磁纬的空间变化项F3项,F3描述了NmF2对地磁纬度的依赖性:
Figure BDA0003453142400000125
Figure BDA0003453142400000126
其中,
Figure BDA0003453142400000127
和λ分别表示观测点的地理纬度和地理经度坐标,
Figure BDA0003453142400000128
Figure BDA0003453142400000129
分别表示北极磁极点的纬度和经度坐标,
Figure BDA00034531424000001210
表示观测点的地磁纬度,c9为模型系数;地磁极点会随时间变化而变化,因此该位置不是固定的常数,可通过国际参考地磁场模型(InternationalGeomagnetic Reference Field model,IGRF)获得;
步骤B6,计算电离层NmF2随磁赤道驼峰区的空间变化项F4,利用高斯函数描述NmF2与磁赤道两侧的驼峰区位置之间的联系:
Figure BDA00034531424000001211
Figure BDA00034531424000001212
其中:
Figure BDA00034531424000001213
Figure BDA00034531424000001214
分别表示北侧和南侧驼峰区磁纬度,分别设置为16°N和15°S,σc表示高斯半峰宽,LT表示地方时,c10—c11为模型系数;
步骤B7,计算电离层NmF2随太阳活动的变化项F5,F5描述了NmF2与太阳射电通量F10.7指数之间的关联:
Figure BDA0003453142400000131
其中,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
Figure BDA0003453142400000132
Figure BDA0003453142400000133
其中,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,计算太阳赤经δ:
Figure BDA0003453142400000141
其中:ζ=180/π,DOY表示年积日;
步骤C3,计算电离层hmF2随地方时LT的时间变化项R1,R1描述了hmF2存在的日、半日和1/3日三个周期变化:
R1=1+ρ1Z1+Z22cosVD3cosVSD4sinVSD5cosVTD6sinVTD7cosVTD)
Figure BDA0003453142400000142
Figure BDA0003453142400000143
Figure BDA0003453142400000144
其中,LT表示地方时,VD、VSD、VTD分别表示三个谐波成分的角相位;Z1和Z2表示hmF2与太阳天顶角χ之间的依赖性,其中
Figure BDA0003453142400000145
表示地理纬度,δ为太阳赤经,ρ1—ρ7为模型系数;
步骤C4:计算电离层hmF2随年积日的时间变化项R2项,R2项描述了hmF2对年和半年的周期变化:
R2=1+c8cos(VA)+c9cos(VSA)
Figure BDA0003453142400000146
Figure BDA0003453142400000147
其中,VA表示年周期变化成分,VSA表示半年周期变化成分,DOY表示年积日,ρ8—ρ9为模型系数;
步骤C5,计算电离层hmF2空间变化项R3,计算方法如下:
Figure BDA0003453142400000151
Figure BDA0003453142400000152
其中,σc表示高斯半峰宽,
Figure BDA0003453142400000153
Figure BDA0003453142400000154
分别设置为40°和20°,LT表示地方时,ρ10—ρ11为模型系数,σLT表示当地时14点对应的高斯半峰宽;
步骤C6,计算电离层hmF2随太阳活动的变化项R4,R4描述了hmF2与太阳射电通量F10.7指数之间的关联:
Figure BDA0003453142400000155
其中,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
Figure BDA0003453142400000156
Figure BDA0003453142400000161
其中,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,构建电离层层析成像线性方程,方程具体表示如下:
Figure BDA0003453142400000171
式中,G是地基GNSS观测的射线数目,B是卫星信标观测的射线数目,N是总网格数,向量dd由地基GNSS和卫星信标接收机的TEC观测数据组成;AA是GNSS和卫星信标无线电信号传播路径在网格中的投影矩阵;X是未知电子密度的分布;e是测量与离散化后引入的误差;
步骤F2,对网格内同一高度层的电子密度加入水平约束,以实现对无射线穿越网格内电子密度的调整。水平约束采用双线性插值算法实现:
Figure BDA0003453142400000172
其中,j表示网格编号,N为总网格数,w表示双线性插值所用的加权系数;
步骤F3,对步骤F1和F2中的线性方程组进行整合,同时忽略离散化误差,层析成像转换为求解下列线性方程组的问题:
d=AX
式中,向量d由dd和0值组合而成,A是投影矩阵AA和加权系数矩阵的组合;
步骤F4,读取步骤D3存储的电子密度数据文件,将其设定为电离层电子密度初值X(0)
步骤F5,利用代数重构算法ART进行电离层三维层析成像,ART算法计算方法如下:
Figure BDA0003453142400000181
其中,i为TEC路径编号;j为网格编号;
Figure BDA0003453142400000182
Figure BDA0003453142400000183
分别为迭代k+1次和k次的第j个网格的电子密度;aij为第i条路径在第j个网格内的投影;||ai||为第i条射线路径的总长;di为第i条射线路径的总电子含量;ξ为松弛因子;
步骤F6,ξk采用与迭代轮次k相关的函数,计算方法如下:
Figure BDA0003453142400000184
步骤F7,设置迭代截止门限tol,当tol<10-8时停止迭代,迭代截止门限计算方法为:
Figure BDA0003453142400000185
其中,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数据是否满足以下条件:
Figure FDA0003453142390000011
Figure FDA0003453142390000012
若满足以上条件,则保存该笔数据,否则认为数据存在较大偏差,剔除该笔数据;
步骤B,电离层F2层峰值电子密度NmF2驱动更新:
步骤B1,利用诺伊斯特雷利兹全球电离层峰值密度模型NPDM进行全球NmF2的变化特征建模,模型表述如下:
NmF2NPDM=F1·F2·F3·F4·F5
其中:F1项描述了电离层NmF2的逐日变化;F2项描述了电离层NmF2的年变化;F3描述了电离层NmF2随磁纬的变化;F4描述了电离层NmF2随磁赤道驼峰区的变化;F5描述了电离层NmF2随太阳活动的变化;
步骤B2,计算太阳赤经δ:
Figure FDA0003453142390000021
其中:ζ=180/π,DOY表示年积日;
步骤B3,计算电离层NmF2随地方时LT的时间变化项F1,F1描述了NmF2存在的日、半日和1/3日三个周期变化:
F1=P1+P2(c1cosVD+c2cosVSD+c3sinVSD+c4cosVTD+c5sinVTD)+c6BO
Figure FDA0003453142390000022
Figure FDA0003453142390000023
Figure FDA0003453142390000024
Figure FDA0003453142390000025
Figure FDA0003453142390000026
其中,LT表示地方时,VD、VSD、VTD分别表示三个不同谐波成分的角相位;P1和P2表示NmF2与太阳天顶角χ之间的依赖性;
Figure FDA0003453142390000027
表示地理纬度,单位弧度;
Figure FDA0003453142390000028
表示地理纬度,单位度;δ为太阳赤经,BO表示中纬区域电离层夏季咬蚀效应;LTBO表示咬蚀效应出现的地方时;c1—c6为模型系数;
步骤B4,计算电离层NmF2随年积日的时间变化项F2项,F2项描述了NmF2对年和半年的周期变化:
F2=1+c7cos(VA)+c8cos(VSA)
Figure FDA0003453142390000029
Figure FDA00034531423900000210
其中,VA表示年周期变化成分,VSA表示半年周期变化成分,c7—c8为模型系数;
步骤B5,计算电离层NmF2随磁纬的空间变化项F3项,F3描述了NmF2对地磁纬度的依赖性:
Figure FDA0003453142390000031
Figure FDA0003453142390000032
其中,
Figure FDA0003453142390000033
和λ分别表示观测点的地理纬度和地理经度坐标,
Figure FDA0003453142390000034
Figure FDA0003453142390000035
分别表示北极磁极点的纬度和经度坐标,
Figure FDA0003453142390000036
表示观测点的地磁纬度,c9为模型系数;地磁极点会随时间变化而变化,位置不是固定的常数,通过国际参考地磁场模型获得;
步骤B6,计算电离层NmF2随磁赤道驼峰区的空间变化项F4,利用高斯函数描述NmF2与磁赤道两侧的驼峰区位置之间的联系:
Figure FDA0003453142390000037
Figure FDA0003453142390000038
其中:
Figure FDA0003453142390000039
Figure FDA00034531423900000310
分别表示北侧和南侧驼峰区磁纬度,分别设置为16°N和15°S,σc表示高斯半峰宽,LT表示地方时,c10—c11为模型系数;
步骤B7,计算电离层NmF2随太阳活动的变化项F5,F5描述了NmF2与太阳射电通量F10.7指数之间的关联:
Figure FDA00034531423900000311
其中,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
Figure FDA0003453142390000041
Figure FDA0003453142390000042
其中,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,计算太阳赤经δ:
Figure FDA0003453142390000043
其中:ζ=180/π,DOY表示年积日;
步骤C3,计算电离层hmF2随地方时LT的时间变化项R1,R1描述了hmF2存在的日、半日和1/3日三个周期变化:
R1=1+ρ1Z1+Z22cosVD3cosVSD4sinVSD5cosVTD6sinVTD7cosVTD)
Figure FDA0003453142390000051
Figure FDA0003453142390000052
Figure FDA0003453142390000053
其中,LT表示地方时,VD、VSD、VTD分别表示三个谐波成分的角相位;Z1和Z2表示hmF2与太阳天顶角χ之间的依赖性,其中
Figure FDA0003453142390000054
表示地理纬度,δ为太阳赤经,ρ1—ρ7为模型系数;
步骤C4:计算电离层hmF2随年积日的时间变化项R2项,R2项描述了hmF2对年和半年的周期变化:
R2=1+c8cos(VA)+c9cos(VSA)
Figure FDA0003453142390000055
Figure FDA0003453142390000056
其中,VA表示年周期变化成分,VSA表示半年周期变化成分,DOY表示年积日,ρ8—ρ9为模型系数;
步骤C5,计算电离层hmF2空间变化项R3,计算方法如下:
Figure FDA0003453142390000057
Figure FDA0003453142390000058
其中,σc表示高斯半峰宽,
Figure FDA0003453142390000059
Figure FDA00034531423900000510
分别设置为40°和20°,LT表示地方时,ρ10—ρ11为模型系数,σLT表示当地时14点对应的高斯半峰宽;
步骤C6,计算电离层hmF2随太阳活动的变化项R4,R4描述了hmF2与太阳射电通量F10.7指数之间的关联:
Figure FDA0003453142390000061
其中,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
Figure FDA0003453142390000062
Figure FDA0003453142390000063
其中,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,构建电离层层析成像线性方程,方程具体表示如下:
Figure FDA0003453142390000071
式中,G是地基GNSS观测的射线数目,B是卫星信标观测的射线数目,N是总网格数,向量dd由地基GNSS和卫星信标接收机的TEC观测数据组成;AA是GNSS和卫星信标无线电信号传播路径在网格中的投影矩阵;X是待求解的电子密度分布;e是测量与离散化后引入的误差;
步骤F2,对网格内同一高度层的电子密度加入水平约束,水平约束采用双线性插值算法实现:
Figure FDA0003453142390000072
其中,j表示网格编号,N为总网格数,w表示双线性插值所用的加权系数;
步骤F3,对步骤F1和F2中的线性方程组进行整合,同时忽略离散化误差,层析成像转换为求解下列线性方程组的问题:
d=AX
式中,向量d由dd和0值组合而成,A是投影矩阵AA和加权系数矩阵的组合;
步骤F4,读取步骤D3存储的电子密度数据文件,将其设定为电离层电子密度初值X(0)
步骤F5,利用代数重构算法ART进行电离层三维层析成像,ART算法计算方法如下:
Figure FDA0003453142390000081
其中,i为TEC路径编号;j为网格编号;
Figure FDA0003453142390000082
Figure FDA0003453142390000083
分别为迭代k+1次和k次的第j个网格的电子密度;aij为第i条路径在第j个网格内的投影;||ai||为第i条射线路径的总长;di为第i条射线路径的总电子含量;ξ为松弛因子;
步骤F6,ξk采用与迭代轮次k相关的函数,计算方法如下:
Figure FDA0003453142390000084
步骤F7,设置迭代截止门限tol,当tol<10-8时停止迭代,迭代截止门限计算方法为:
Figure FDA0003453142390000085
其中,N是总网格数;
步骤F8,ART算法迭代终止后,输出最终的X(k+1)到文件中并存储,即为最终的电离层电子密度层析成像结果。
CN202111678437.7A 2021-12-31 2021-12-31 一种基于多源数据驱动的电离层层析成像方法 Active CN114384564B (zh)

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)

* Cited by examiner, † Cited by third party
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)

* Cited by examiner, † Cited by third party
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

Patent Citations (3)

* Cited by examiner, † Cited by third party
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)

* Cited by examiner, † Cited by third party
Title
欧明;甄卫民;徐继生;刘钝;於晓;: "地基GPS与掩星联合的电离层层析成像方法研究", 全球定位系统, no. 05, 31 October 2014 (2014-10-31) *

Cited By (3)

* Cited by examiner, † Cited by third party
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