三频信标机数据处理方法及系统与流程

1.本发明涉及卫星遥感技术领域，特别是涉及一种三频信标机数据处理方法及系统。


背景技术：

2.““
张衡一号”卫星是中国第一颗自主研制的地球物理场卫星，卫星上有八大载荷。“张衡一号”卫星于2018年2月2日成功发射，该卫星搭载的全球导航卫星系统无线电掩星(gro)接收机具备了全球定位系统(gps)卫星和北斗导航系统(bd)卫星的掩星观测功能。


技术实现要素：

3.本发明要解决的技术问题是提供一种三频信标机数据处理方法及系统，能够将三频信标机的原始探测数据进行处理，生成实际可用的数据产品。
4.为解决上述技术问题，本发明提供了一种三频信标机数据处理方法，所述方法包括：地面三频信标接收机跟踪锁定星载三频信标信号后，输出三个频段信号的正交分量q和同向分量i观测数据和信噪比数据；对三频信标接收机的i/q观测数据分别进行差分相位计算和信号强度计算，可计算得到vhf&uhf、uhf&l的差分相位值和三个频段各自的信号强度值；对计算得到的差分相位值进行相位连接处理，每秒计算出一组相对电离层总电子含量tec值，对计算得到的信号强度值进行幅度闪烁指数计算，每1秒计算出一组闪烁指数，进行时空匹配，计算得到各组数据在观测时间所对应的空间位置；通过对多个台站的观测数据进行集中处理，利用多站法每秒钟计算一个电离层绝对tec数据，再利用层析成像算法对单个链路进程反演得到二维电子密度剖面数据，并提取nmf2异常。
5.在一些实施方式中，三个频段信号的正交分量q和同向分量i观测数据和信噪比数据被称为0级数据；vhf&uhf、uhf&l的差分相位值和三个频段各自的信号强度值被称为1级数据；对计算得到的差分相位值进行相位连接处理，计算出的相对电离层总电子含量值(相对tec)，和根据信号强度值进行幅度闪烁指数计算得到的闪烁指数被称为2级数据；利用层析成像算法对单个链路进程反演得到的二维电子密度剖面数据，及提取到的nmf异常被称为3级数据。
6.在一些实施方式中，0级数据包括：地面站接收机产出的带有时间信息的原始观测数据；1级数据包括：对各单站0级数据进行正交解算后得到的按时间排列的相位数据和幅度数据；2级数据包括：对各单站1级数据进行相位连接和闪烁指数计算的带有轨道信息的数据；3级数据包括：在2级数据的基础上，对多站数据进行反演处理生成一条观测链路的绝对tec和二维剖面产品。
7.在一些实施方式中，0级数据包括：vhf、uhf、l三个频段的i/q观测数据、信噪比数据文件；1级数据包括：科学数据、图像产品和数据处理报告；2级数据包括：科学数据、图像产品和数据处理报告；3级数据包括：科学数据及其图像和数据处理报告。
8.在一些实施方式中，对三频信标接收机的i/q观测数据分别进行差分相位计算和
信号强度计算，可计算得到vhf&uhf、uhf&l的差分相位值和三个频段各自的信号强度值，包括：差分相位计算，以及信号强度计算。
9.在一些实施方式中，差分相位计算，包括：读取0级vhf频段i路、q路数据；计算vhf/uhf差分相位；读取0级l频段i路、q路数据；计算l/uhf差分相位。
10.在一些实施方式中，信号强度计算，包括：读取0级vhf、uhf、l频段的i路、q路数据；计算vhf、uhf、l频段的信号强度。
11.在一些实施方式中，对计算得到的差分相位值进行相位连接处理，每秒计算出一组相对电离层总电子含量tec值，对计算得到的信号强度值进行幅度闪烁指数计算，每1秒计算出一组闪烁指数，进行时空匹配，计算得到各组数据在观测时间所对应的空间位置，包括：相对tec计算，以及闪烁指数计算。
12.在一些实施方式中，通过对多个台站的观测数据进行集中处理，利用多站法每秒钟计算一个电离层绝对tec数据，再利用层析成像算法对单个链路进程反演得到二维电子密度剖面数据，并提取nmf2异常，包括：获取本站链对应的所有接收站的相对tec数据，利用多站法，计算出每个接收站的积分常数，将积分常数加至接收机得到的的相对tec，得到各个接收站的绝对tec；将站链上空100km～500km范围按照高度间隔20km，纬度间隔0.5
°
进行间隔划分，通过反演过程计算划分区域的电子密度剖面。
13.此外，本发明还提供了一种三频信标机数据处理系统，所述系统包括：一个或多个处理器；存储装置，用于存储一个或多个程序，当所述一个或多个程序被所述一个或多个处理器执行，使得所述一个或多个处理器实现根据前文所述的三频信标机数据处理方法。
14.本发明提供的三频信标机数据处理方法及系统，通过对原始探测数据的差分相位计算和信号强度计算等一系列运算，生成数据精度及格式均符合要求的数据产品。
附图说明
15.上述仅是本发明技术方案的概述，为了能够更清楚了解本发明的技术手段，以下结合附图与具体实施方式对本发明作进一步的详细说明。
16.图1是三频信标机系统组成框图；
17.图2是三频信标观测系统；
18.图3是0级数据命名规则的原理示意图；
19.图4是1
‑
4级数据命名规则的原理示意图；
20.图5是科学数据处理报告命名的原理示意图；
21.图6是科学数据产品图像命名的原理示意图；
22.图7是三频信标机1级图像产品示意图；
23.图8是三频信标机2级图像产品示意图；
24.图9是三频信标3级数据示意图；
25.图10是三频信标3级nmf2动态变化图；
26.图11是三频信标机0～3级数据处理流程图；
27.图12是三频信标0级
‑
1级数据处理流程图；
28.图13是三频信标1级
‑
2级数据处理流程图；
29.图14是相对tec计算处理流程；
30.图15是闪烁指数计算流程；
31.图16是2级
‑
3级数据处理流程；
32.图17是求积分常数计算示意图；
33.图18是k2计算流程图；
34.图19是电子密度反演网格示意图。
具体实施方式
35.以下结合附图对本发明的优选实施例进行说明，应当理解，此处所描述的优选实施例仅用于说明和解释本发明，并不用于限定本发明。
36.首先介绍下文中需要使用的缩略语。本技术下文中使用的缩略语如表1所示。
37.表1缩略语
[0038][0039]
1三频信标机简介
[0040]
1.1系统工作组成与框图
[0041]
三频信标观测系统利用卫星平台能够对电离层进行大范围、高时间分辨率、高水平分辨率的扫描，并能对f2层以上电离层的状态以及对小尺度的电离层不均匀体结构进行测量，能够获取高精度电离层tec及电离层电子密度等信息。
[0042]
三频信标观测系统由星载发射系统和地面接收系统组成，如图1所示。其中，星载发射系统由三品信标发射机、三频信标发射天线以及它们之间的射频电缆组成；地面接收系统由接收天线、接收主机、gps卫星授时、显示控制四个部分组成。
[0043]
1.2工作原理
[0044]
三频信标系统观测原理如图2所示。通过对所要探测电离层区域的全方位断层扫描，进而利用扫描投影数据重建该区域内电离层电子密度等参量的空间分布。三频信标观测系统通过星载发射系统向地面发射三个频率的相干信号，三频信标地面接收系统接收上述信号，并进行相关分析计算，反演出相应位置的电子密度分布，从而实现对电离层的断层扫描。
[0045]
星载发射机内部的恒温晶振为发射机提供高稳定的参考信号。发射机运行以后，经过频率源模块进行频率变换，产生vhf频段(中心频率：150.012mhz)、uhf频段(中心频率：
400.032mhz)及l频段(中心频率：1066.752mhz)的相干射频信号。这三个射频信号分别经过正交电桥，分路产生等幅的6路正交信号，经功率放大后，分别接入发射天线的i路输入端和q路输入端，向空间进行辐射圆极化波，其穿越电离层后信号的接收任务由地面的三频信标测量站来完成。地面站链上的各台接收机通过天线接收信号，经过三波段信标信号相干处理等运算后，产出原始观测数据。
[0046]
三频信标接收机系统主要由接收天线、接收主机、gps卫星授时、显示控制四部分构成。接收天线主要完成对卫星信标信号的接收和初步预处理。接收天线使用的是正交馈电的正交偶极子，并采用圆极化的方式接收信号。天线处理单元主要完成信标信号的滤波、合路以及功率放大功能。接收机主要完成信标信号的核心处理过程，主要包括卫星信标信号的捕获、卫星信标信号的高精度跟踪、三波段信标信号的相干处理等内容。接收机主要由射频前端和中频信号处理两部分组成，其中射频前端单元完成信标信号的相干下变频变换，中频信号处理单元完成核心算法的实现。gps卫星授时部分完成对接收机的精确授时，并可实现多台接收机之间的时间同步功能。
[0047]
1.3数据产出
[0048]
三频信标机观测数据由地面站接收机产出。接收机通过星历文件计算出卫星飞越其上空的时间，提前开机准备接收卫星信号。卫星从地面站上空飞过时，每轨可由接收机产出一个十进制的i/q观测数据文件，以及一个十进制的信噪比数据文件。
[0049]
2标准数据产品
[0050]
2.1分级定义
[0051]
按照《卫星对地观测数据产品分类分级规则》(gb/t 32453
‑
2015)，以及《电磁监测卫星数据产品分级及定义》，三频信标观测数据包括了0级，1级、2级和3级数据产品，三频信标标准数据产品定义见表2。
[0052]
表2 tbb数据产品分级定义及数据产品一览表
2payloadid仪器代码 3orbitnum轨道号 4orbittype轨道类别1：升轨，0：降轨
[0068]
表4三频信标机1级数据说明表
[0069][0070][0071]
(2)1级图像产品
[0072]
三频信标接收机的1级图像产品记录每条轨道的vhf、uhf、l频段的差分相位、信号幅度随时间的变化趋势，共包含6组数据曲线：
[0073]
◆
vhf&uhf差分相位
‑
t
[0074]
◆
uhf&uhf差分相位
‑
t
[0075]
◆
l&uhf差分相位
‑
t
[0076]
◆
vhf信号幅值
‑
t
[0077]
◆
uhf信号幅值
‑
t
[0078]
◆
l信号幅值
‑
t。
[0079]
三频信标接收机1级图像产品示意图如图7所示。
[0080]
(3)1级数据处理报告
[0081]
三频信标机1级数据处理报告组成要素：
[0082]
(1)处理版本号：v0.1
[0083]
(2)起始时间:yyyymmdd hh:mm:ss.zzz
[0084]
(3)输入数据：0级数据文件名
[0085]
(4)辅助数据
[0086]
①
功放温度：正常/不正常；标注出不正常时间段；
[0087]
②
12v、9.5v、
±
5v电压：正常/不正常；标注出不正常时间段；
[0088]
(5)处理过程
[0089]
①
正常
[0090]
②
缺数
[0091]
③
数据文件损坏
[0092]
(6)结束时间：yyyymmdd hh:mm:ss.zzz
[0093]
(7)输出产品：1级数据产品文件名
[0094]
2.2.3 2级数据产品
[0095]
2级数据是对1级数据vhf、uhf、l频段的信号相位数据进行相位连接等数学处理而求得的相对tec数据，以及利用vhf、uhf、l频段的信号幅度数据计算出来的电离层闪烁指数数据。2级数据产品包括：科学数据、图像产品和数据处理报告。
[0096]
(1)科学数据
[0097]
三频信标机2级科学数据文件包含以下内容：
[0098]
◆
vhf&uhf频段相对tec值
[0099]
◆
l&uhf频段相对tec值
[0100]
◆
vhf频段闪烁指数值
[0101]
◆
uhf频段闪烁指数值
[0102]
◆
l频段闪烁指数值
[0103]
三频信标机2级数据格式见表5、表6所示。
[0104]
表5三频信标机2级数据文件属性说明
[0105]
序号属性名称属性内容备注1softversion程序版本号 2payloadid仪器代码 3orbitnum轨道号 4orbittype轨道类别1：升轨，0：降轨
[0106]
表6三频信标机2级数据表格式说明
[0107][0108]
(2)2级图像产品
[0109]
三频信标机2级图像产品包含以下5组曲线：
[0110]
◆
vhf&uhf频段相对tec值
‑
t
[0111]
◆
l&uhf频段相对tec值
‑
t
[0112]
◆
vhf频段闪烁指数
‑
t
[0113]
◆
uhf频段闪烁指数
‑
t
[0114]
◆
l频段闪烁指数
‑
t
[0115]
三频信标机2级图像产品示意图如图8所示。
[0116]
(3)2级数据处理报告
[0117]
三频信标机2级数据处理报告文件组成要素：
[0118]
(1)处理软件版本号：v0.1
[0119]
(2)起始时间：yyyymmdd hh:mm:ss.zzz
[0120]
(3)输入数据：1级数据文件名
[0121]
(4)辅助数据
[0122]
(5)处理过程
[0123]
①
处理过程正常
[0124]
②
异常现象及可能的原因
[0125]
(6)结束时间：yyyymmdd hh:mm:ss.zzz
[0126]
(7)输出数据：2级数据文件名
[0127]
2.2.4 3级数据产品
[0128]
3级数据是在2级数据的基础上，进行反演处理生成一条观测链路的绝对tec和二维剖面产品。根据电子密度二维分布结果，提取nme、nmf1、nmf2异常。3级数据产品包括：科学数据及其图像和数据处理报告。
[0129]
科学数据
[0130]
3级数据产品格式详见表7、表8。
[0131]
表7三频信标3级数据文件属性说明
[0132]
序号属性名称属性内容备注1payloadid仪器代码 2orbitnum轨道号 3orbittype轨道类别升轨、降轨4softversion程序版本号vr0.1
[0133]
表8三频信标3级数据表格式说明
[0134][0135][0136]
(2)图像产品
[0137]
3级数据产品示意图如图9、图10所示。
[0138]
(2)3级数据处理报告
[0139]
3级数据处理报告文件类型为txt文件。组成要素：
[0140]
(1)处理软件版本号：v0.1
[0141]
(2)起始时间：yyyymmdd hh:mm:ss.zzz
[0142]
(3)输入数据：(a)2级数据文件名；(b)前5条重访轨道的3级数据，列出其文件名
[0143]
(4)辅助数据
[0144]
①
地震事件：该轨道时间之前1天内、100km范围6级以上地震记录
[0145]
②
磁情指数：该轨对应时间的kp、dst等
[0146]
③
iri模型
[0147]
(5)处理过程
[0148]
①
处理过程正常
[0149]
②
异常现象及可能的原因
[0150]
(6)结束时间：yyyymmdd hh:mm:ss.zzz
[0151]
(7)输出数据：3级数据文件名
[0152]
3数据处理流程
[0153]
地面三频信标接收机跟踪锁定星载三频信标信号后，输出三个频段信号的正交分量q和同向分量i观测数据和信噪比数据(level 0数据)。对三频信标接收机的i/q观测数据分别进行差分相位计算和信号强度计算，可计算得到vhf&uhf、uhf&l的差分相位值和三个频段各自的信号强度值(level 1数据)。对计算得到的差分相位值进行相位连接处理，每秒计算出一组相对电离层总电子含量tec值；对计算得到的信号强度值进行幅度闪烁指数计算，每1秒计算出一组闪烁指数。进行时空匹配，计算得到各组数据在观测时间所对应的空间位置(level 2数据)；通过对多个台站的观测数据进行集中处理，利用多站法每秒钟计算一个电离层绝对tec数据，再利用层析成像算法对单个链路进程反演得到二维电子密度剖面数据，并提取nmf2异常(level 3)。三频信标系统各级数据处理流程图见图11所示。
[0154]
4数据处理方法
[0155]
4.1 0级数据生成1级数据
[0156]
三频信标0级数据生成1级数据的处理流程如图12所示。
[0157]
4.1.1差分相位计算
[0158]
第1步：读取0级vhf频段i路、q路数据
[0159]
三频信标0级数据主要是vhf与uhf、l与uhf频段两两组合，进行相位差分处理后得到的i路(同向)和q路(正交)两组数据。其中uhf是接收机锁定信号时进行差分处理的基准频率，因此得到的有效差分数据只有vhf/uhf和l/uhf两种组合。每秒钟50组数据。
[0160]
第2步：计算vhf/uhf差分相位
[0161]
按照公式(1)进行反正切处理，计算得到vhf/uhf频段对应的差分相位。
[0162]
φ(t)＝arctan(q/i)
ꢀꢀꢀ
(1)
[0163]
第3步：读取0级l频段i路、q路数据
[0164]
读取l频段i路和q路数据，每秒钟50组数据。
[0165]
第4步：计算l/uhf差分相位
[0166]
按照公式(1)进行反正切计算，得到l/uhf频段对应的差分相位数据。
[0167]
4.1.2信号强度计算
[0168]
第1步：读取0级vhf、uhf、l频段的i路、q路数据
[0169]
分别读取三频信标0级vhf、uhf、l频段的i路、q路数据。
[0170]
第2步：计算vhf、uhf、l频段的信号强度
[0171]
按照公式(2)进行信号强度计算。计算结果为功率值，分别得到vhf、uhf、l频段对应的三组信号强值数据，单位为dbm。
[0172][0173]
其中，10lg(i2 q2)的处理是把接收机处理产出的相对电压值转换为功率值，231(db)是通道增益量(包含天线增益、射频增益和算法增益)。
[0174]
4.2 1级数据生成2级数据
[0175]
三频信标1级数据生成2级数据的处理流程如图13所示。
[0176]
4.2.1相对tec计算
[0177]
tec是指电离层电子密度沿卫星
‑
接收机链路的路径积分量。相干信标接收机所记录的差分多普勒相位限制在0～360
°
之间，因此必须将数据连接起来，才能变成相应的相对tec。相位连接方法是在数据相位0～360
°
中取一个合适的门限值对相位进行判断和数据连接，具体实现步骤如图14所示。
[0178]
具体处理方法流程为：
[0179]
第1步：读取差分相位数据
[0180]
分别读取vhf&uhf、uhf&l频段的差分相位数据。
[0181]
第2步：取初始值，设置相位门限
[0182]
取t＝1，k(1)＝1，d＝300。
[0183]
第3步：差分计算
[0184]
对于有n个点的相位数据，按照公式(3)所示对相位数据进行差分计算。
[0185]
δ(t)＝φ(t 1)
‑
φ(t)，t＝1,2,....,n
‑1ꢀꢀꢀ
(3)
[0186]
第4步：相位翻转个数计算
[0187]
将各数据点的δ值与相位门限值d按照公式(4)的关系进行比较，得到各相位数据在连接时所取2π的翻转个数k。
[0188][0189]
第5步：相位连接处理
[0190]
按照公式(5)依次对各数据点的差分相位序列进行相位值的增加和减少，即相位连接。
[0191]
φ'
c
(t)＝φ(t) 2kπ
ꢀꢀꢀ
(5)
[0192]
第6步：计算最小相位值
[0193]
在连接完成后，按公式(6)计算相位序列中的最小相位值z。
[0194]
z＝min(φ'
c
(t))
ꢀꢀꢀ
(6)
[0195]
第7步：最小相位归零计算
[0196]
φ
c
(t)＝φ'
c
(t)
‑
z
ꢀꢀꢀ
(7)
[0197]
其中φ'
c
(t)为连接起来的相位序列，按照公式(7)依次减掉最小相位值z，最终得到连接完成的差分多普勒相位序列φ
c
(t)。
[0198]
第8步：计算相关系数c
d
[0199][0200]
按照公式(8)计算相关系数c
d
。其中，c＝299792458m/s，f
r
＝16.668mhz，m
i
表示倍频数，如果取vhf/uhf组合，则取m1＝9，m2＝24；如果取uhf/l组合，则取m1＝24，m2＝64。
[0201]
第9步：计算相对tec值
[0202]
按照公式(4
‑
9)可计算得到电离层相对tec值。tec值的单位为tecu(1tecu＝1016el.m
‑
2)。
[0203]
tec
r
＝c
d
×
φ
c
(t)
ꢀꢀꢀ
(9)
[0204]
4.2.2闪烁指数计算
[0205]
三频信标机幅度闪烁指数s4的过程如图15所示。
[0206]
第1步：读取信号强度值
[0207]
分别读取vhf、uhf、l频段的信号强度值数据。
[0208]
第2步：数据分组
[0209]
以每50个数(1秒钟长度)为1组进行分组。
[0210]
第3步：计算各组数据平方的均值
[0211][0212]
按照公式(10)计算得到每组数据各个数值的平方的平均值。
[0213]
第4步：计算各组数据平均值的平方
[0214][0215]
按照公式(11)计算得到每组数据的平均值的平方。
[0216]
第5步：计算闪烁指数值
[0217][0218]
按照公式(12)每1秒钟计算得到一个闪烁指数值。一般0.1≤s4<0.3被认为是电离层弱闪烁，0.3≤s4<0.6表示电离层中等闪烁，s4>0.6表示电离层强闪烁。
[0219]
4.3 2级数据生成3级数据
[0220]
三频信标2级数据生成3级数据的处理流程如图16所示。
[0221]
(1)获取本站链对应的所有接收站的相对tec数据，利用多站法，计算出每个接收站的积分常数，将积分常数加至接收机得到的的相对tec，即得到各个接收站的绝对tec。
[0222]
绝对tec具体计算步骤如下：
[0223]
首先利用双频信标计算绝对tec。
[0224]
如图17，设接收机a到卫星s的斜tec为i
s1
，其投影的垂向tec为i
v1
，则二者具有如下关系：
[0225][0226]
其中，d1称为路径s1上斜向全电子含量与垂直全电子含量的转换函数，χ
a
为路径s1与接收站a垂线方向的夹角。卫星过境期间斜向全电子含量i
s1
与差分多普勒相位关系式表示为：
[0227][0228]
其中i
a
(t
k
)为接收机a在t
k
(k＝1，
…
，n)时刻观测给出的差分相对tec；φ1(0)为待求接收站a的未知积分常数。其值等于1.3
×
10,3，其中c为光速，f0＝16.668mhz，n1＝3,n2＝8；φ
d1
(t
k
)是接收机a在各时刻观测对应的差分多普勒相位；公式(14)可表示为
[0229][0230]
其中，d
1k
和i
v1k
分别为第k个观测点与a对应射线的转换函数和垂直tec。同理，对于接收站b有：
[0231][0232]
公式(15)、(16)可以变为
[0233][0234][0235]
近似地认为，对于地面不同位置的两个接收站，两个接收站的垂向tec相等。因此对于双站法i
v1k
≈i
v2k
，利用最小平方误差限定i
v1k
和i
v2k
两者间的误差，双站法的目标函数可表示下式：
[0236][0237]
通过对e求φ1(0)和φ2(0)的偏导数，再令其偏微分值为零，配合行列式求解。即联立方程立方程即可求得两站的未知积分常数φ1(0)和φ2(0)。
[0238]
从而带入公式(14)可由观测得到的差分相对tec得到各个站的双频绝对tec如下式所示：
[0239][0240][0241]
推广到m个接收站，也可按照相同的概念求解多个站的初始相位，即所谓的多站法(multi
‑
station method)。多站法的目标函数为公式(4
‑
20)。
[0242][0243]
通过求解e的极小值也可得到各个站的未知相位积分常数φ
i
(0)(i＝1,
···
,m)，从而求得各个站的双频绝对tec。
[0244][0245]
②
利用双频辅助三频求解三频信标绝对tec
[0246]
三频信标与双频信标相比多了一个载波两个频率的载波可做如下差分：
[0247][0248]
在三频信标中有:
[0249]
f1＝9f0＝150.012mhz
[0250]
f2＝24f0＝400.032mhz
[0251]
f3＝64f0＝1066.752mhz
[0252]
因此三组载波利用差分多普勒技术可以得到三个差分相位表达式(23)：
[0253][0254]
式中，f0＝16.668mhz，n1＝9,n2＝24,n3＝64，δp
12
、δp
13
、δp
23
为差分相位的绝对值，包含接收机给出的差分多普勒相位和相位积分常数。
[0255]
则：
[0256][0257]
由公式(24)可得：
[0258]
δp
12
(n
12
n
22
)＝δp
13
n
22
ꢀꢀꢀ
(25)
[0259]
上式可改写为：
[0260]
(k
12
δφ
12
)(n
12
n
22
)＝(k
13
δφ
13
)n
22
或
[0261]
k
12
(n
12
n
22
)
‑
k
13
n
22
＝δφ
13
n
22
‑
δφ
12
(n
12
n
22
)
ꢀꢀꢀ
(26)
[0262]
其中，k
12
和k
13
是δp
12
和δp
13
整数部分，δφ
12
和δφ
13
是其小数部分，由接收机数据给出。基于素数原理，有：
[0263]
x
12
(n
12
n
22
)
‑
x
13
n
22
＝1
ꢀꢀꢀ
(27)
[0264]
其中，n1和n2互质，联立公式(26)、(27)可解出k
12
和k
13
，其表达式为：
[0265]
k
12
＝[δφ
13
n
22
‑
δφ
12
(n
12
n
22
)]x
12
kn
22
[0266]
k
13
＝[δφ
13
n
22
‑
δφ
12
(n
12
n
22
)]x
13
k(n
12
n
22
)
ꢀꢀꢀ
(28)
[0267]
其中k为一整数，则由k
12
和k
13
可求得δp
12
和δp
13
，再把δp
13
或δp
12
带入式(23)并整理可得到：
[0268][0269]
在三频信标中n1＝3，n2＝8，f0＝16.668mhz。带入上式得到
[0270]
tec＝8.3165
×
10
16
×
[(δφ
13
x
12
‑
δφ
12
x
13
)
mod1
k2]
ꢀꢀꢀ
(30)
[0271]
其中8.3165
×
1016(单位：el/m2)称为三频信标测量tec时的模糊系数，k2为一正整数，它与双频信标的相位积分常数一样，为待定值。
[0272]
由公式(19)、(30)，我们可以将双频和三频tec的求解简化表示为：
[0273]
tec＝k1×
1.2995
×
10
15
φ1ꢀꢀꢀ
(31)
[0274]
tec＝k2×
8.3165
×
10
16
φ2ꢀꢀꢀ
(32)
[0275]
其中，k1和k2分别表示双频和三频相位积分常数，其值均为整数，φ1和φ2分别表示双频和三频差分相位的小数部分。公式(31)通过多站法可以计算出包含误差的相位积分常数k1，再结合地面接收机测得的差分相位小数部分φ1，可以计算出包含误差的绝对tec。由于同一反演区域同一时刻的tec是相同的，将公式(31)计算得到的绝对tec带入公式(32)，即得到k2的计算公式：
[0276][0277]
因为公式(31)计算得到的绝对tec存在误差，所以公式(33)计算得到的k2一般也存在误差，我们需要用三频相位测量结果对其进行精确修正。研究表明，双频测量tec的误差在1～3tecu量级(1tecu＝1.0
×
1016el/m2)，而三频相位模糊系数为8.3tecu(即以8.3tecu为单位分组)，结合k2为整数的基本前提，对k2进行精确求解的算法流程图如图18。
[0278]
(2)将站链上空100km～500km范围按照高度间隔20km，纬度间隔0.5
°
进行间隔划分(图19)，在利用卫星信号观测进行的电离层层析成像过程中，所获得的电离层电子总含量tec是电离层电子密度沿卫星至接收机射线路径上的积分，可表示为：
[0279]
tec＝∫
p
n
e
ds
ꢀꢀꢀ
(34)
[0280]
其中ne表示电子密度分布,p是接收机和卫星间的路径。
[0281]
电子密度反演计算过程如下：
[0282]
①
先利用球谐分析和经验正交函数组合的函数基模型法进行电离层ct。球谐函数和经验正交函数选择用于电子密度的表示：
[0283][0284]
其中θ和分别代表纬度和经度,h代表高度.是归一化的勒让德函数,m,n,k是球谐函数和经验正交函数(eof)的阶数，其中eof通过经验模式或测量的电子密度数据通过奇异值分解得到。
[0285]
根据卫星轨道和接收机位置确定反演区域后，将待反演区域内连续分布的电子密度离散化成n个不同网格，并假设每一个网格内电子密度为恒定不变的。然后确定反演区域内的有效射线，有效射线是指与反演划定区域上、下高度边界(与地球表面同心的球壳)的交点均落在反演区域内的射线。每一个网格内射线穿过的路径投影长度可按如下方法计算：首先将卫星轨道及地面接收站的球坐标转换成为直角坐标(x,y,z)，其转换公式为：
[0286][0287]
其中r＝r
e
h，r
e
表示地球半径。然后利用直线方程将每一条射线通过各个网格上的交点在直角坐标系下的坐标求出。地面接收站与卫星间射线的直线方程可表示为：
[0288][0289]
其中(x
sat
,y
sat
,z
sat
)与(x
sta
,y
sta
,z
sta
)分别表示卫星与三频信标接收机(站)的位置坐标。(x,y,z)表示直线上的任何一点，k表示比例常数，对于每个固定的点(x,y,z)来说，k都是一个可以确定的常量。由上式我们可以得到：
[0290][0291]
若交点位于距地球中心点距离为r的面上，则r已知，因此有：
[0292]
x2 y2 z2＝r2ꢀꢀꢀ
(39)
[0293]
将上式代入其中，再由一元二次方程ak2 bk c＝0求解得到k，代回上式中即可得到相交点的位置坐标，然后再计算
[0294][0295]
求得各网格间的投影长度δl。所有射线在反演区域内各网格间的投影长度的集合即可表示为电离层电子密度反演的投影矩阵。确定投影矩阵以后，通过转换，电离层电子密度反演的问题可以简化为：
[0296]
y＝hx
ꢀꢀꢀ
(41)
[0297]
其中x是基函数的权重系数，绝对tec可以用一列的y来表示.h是由经验正交函数和球谐函数及几何路径投影长度组成的算子。由于卫星信标多站法估计绝对tec存在较大误差，因此我们对以上矩阵进行了修正，我们采用对应的一个接收机tec相互差分的方法来多站法误差的影响，表示为:
[0298]
b＝y
‑
y0＝(h
‑
h0)x＝ax
ꢀꢀꢀ
(42)
[0299]
其中b是tec数据的差分结果，y0是给定接收机的一个参考测量值.其中前向算子a包含了h和参考算子h0。采用截断奇异值分解(tsvd)正则化的方法求解矩阵方程，根据矩阵的奇异值分解理论，可将矩阵a分解为
[0300][0301][0302]
式中：u＝(u1,u2,...,u
n
)，u
i
为左奇异值向量；v＝(v1,v2,...,v
n
)，v
i
为右奇异值向量；σ＝(σ1,σ2,...,σ
n
)，σ
i
为奇异值，σ1≥σ2≥
…
≥σ
n
≥0,k为截断系数。通过式(44)求解得到基函数权重x后，重新带入式(41)便可得到电子密度的初始分布值。
[0303]
②
将函数基模型的反演结果作为迭代初值，采用像素基模型的art算法，对式(42)进行第二次迭代重构，以获取更小尺度的电离层扰动结构，计算方法如下：
[0304][0305]
式中，q表示迭代轮次，具体次序可按顺序或通过随机数的方法产生，i＝mod(q，l) 1,l为矩阵y的行数，t
i
表示矩阵y的第i行元素，r
i
表示r的第i个元素，λ为松弛因子。art算法收敛较快，一般取迭代10～20轮左右即可，也可以在具体执行过程中设置适当的迭代终止阀值。通过以上两步的反演即可获得最终的电离层电子密度分布。
[0306]
以上所述，仅是本发明的较佳实施例而已，并非对本发明作任何形式上的限制，本领域技术人员利用上述揭示的技术内容做出些许简单修改、等同变化或修饰，均落在本发明的保护范围内。