具体实施方式

以下结合附图对本发明的优选实施例进行说明，应当理解，此处所描述的优选实施例仅用于说明和解释本发明，并不用于限定本发明。

首先介绍下文中需要使用的缩略语。本申请下文中使用的缩略语如表1所示。

表1缩略语

1三频信标机简介

1.1系统工作组成与框图

三频信标观测系统利用卫星平台能够对电离层进行大范围、高时间分辨率、高水平分辨率的扫描，并能对F2层以上电离层的状态以及对小尺度的电离层不均匀体结构进行测量，能够获取高精度电离层TEC及电离层电子密度等信息。

三频信标观测系统由星载发射系统和地面接收系统组成，如图1所示。其中，星载发射系统由三品信标发射机、三频信标发射天线以及它们之间的射频电缆组成；地面接收系统由接收天线、接收主机、GPS卫星授时、显示控制四个部分组成。

1.2工作原理

三频信标系统观测原理如图2所示。通过对所要探测电离层区域的全方位断层扫描，进而利用扫描投影数据重建该区域内电离层电子密度等参量的空间分布。三频信标观测系统通过星载发射系统向地面发射三个频率的相干信号，三频信标地面接收系统接收上述信号，并进行相关分析计算，反演出相应位置的电子密度分布，从而实现对电离层的断层扫描。

星载发射机内部的恒温晶振为发射机提供高稳定的参考信号。发射机运行以后，经过频率源模块进行频率变换，产生VHF频段(中心频率：150.012MHz)、UHF频段(中心频率：400.032MHz)及L频段(中心频率：1066.752MHz)的相干射频信号。这三个射频信号分别经过正交电桥，分路产生等幅的6路正交信号，经功率放大后，分别接入发射天线的I路输入端和Q路输入端，向空间进行辐射圆极化波，其穿越电离层后信号的接收任务由地面的三频信标测量站来完成。地面站链上的各台接收机通过天线接收信号，经过三波段信标信号相干处理等运算后，产出原始观测数据。

三频信标接收机系统主要由接收天线、接收主机、GPS卫星授时、显示控制四部分构成。接收天线主要完成对卫星信标信号的接收和初步预处理。接收天线使用的是正交馈电的正交偶极子，并采用圆极化的方式接收信号。天线处理单元主要完成信标信号的滤波、合路以及功率放大功能。接收机主要完成信标信号的核心处理过程，主要包括卫星信标信号的捕获、卫星信标信号的高精度跟踪、三波段信标信号的相干处理等内容。接收机主要由射频前端和中频信号处理两部分组成，其中射频前端单元完成信标信号的相干下变频变换，中频信号处理单元完成核心算法的实现。GPS卫星授时部分完成对接收机的精确授时，并可实现多台接收机之间的时间同步功能。

1.3数据产出

三频信标机观测数据由地面站接收机产出。接收机通过星历文件计算出卫星飞越其上空的时间，提前开机准备接收卫星信号。卫星从地面站上空飞过时，每轨可由接收机产出一个十进制的I/Q观测数据文件，以及一个十进制的信噪比数据文件。

2标准数据产品

2.1分级定义

按照《卫星对地观测数据产品分类分级规则》(GB/T 32453-2015)，以及《电磁监测卫星数据产品分级及定义》，三频信标观测数据包括了0级，1级、2级和3级数据产品，三频信标标准数据产品定义见表2。

表2 TBB数据产品分级定义及数据产品一览表

2.2各级数据产品介绍

2.2.2 1级数据产品

1级数据是对0级数据VHF、UHF、L三个频点的I/Q数据进行计算后得到的单站相位数据和幅度数据。1级数据产品包括：科学数据、图像产品和数据处理报告。

(1)科学数据

三频信标机1级科学数据文件包含以下物理量数据：

◆VHF&UHF频段的差分相位数据

◆UHF&UHF频段的差分相位数据

◆L&UHF频段的差分相位数据

◆VHF频段的信号幅度数据

◆UHF频段的信号幅度数据

◆L频段的信号幅度数据

三频信标机1级数据文件格式见表3、表4所示。

表3三频信标机1级文件属性说明表

表4三频信标机1级数据说明表

(2)1级图像产品

三频信标接收机的1级图像产品记录每条轨道的VHF、UHF、L频段的差分相位、信号幅度随时间的变化趋势，共包含6组数据曲线：

◆VHF&UHF差分相位-t

◆UHF&UHF差分相位-t

◆L&UHF差分相位-t

◆VHF信号幅值-t

◆UHF信号幅值-t

◆L信号幅值-t。

三频信标接收机1级图像产品示意图如图7所示。

(3)1级数据处理报告

三频信标机1级数据处理报告组成要素：

(1)处理版本号：V0.1

(2)起始时间:yyyymmdd HH:MM:SS.ZZZ

(3)输入数据：0级数据文件名

(4)辅助数据

①功放温度：正常/不正常；标注出不正常时间段；

②12V、9.5V、±5V电压：正常/不正常；标注出不正常时间段；

(5)处理过程

①正常

②缺数

③数据文件损坏

(6)结束时间：yyyymmdd HH:MM:SS.ZZZ

(7)输出产品：1级数据产品文件名

2.2.3 2级数据产品

2级数据是对1级数据VHF、UHF、L频段的信号相位数据进行相位连接等数学处理而求得的相对TEC数据，以及利用VHF、UHF、L频段的信号幅度数据计算出来的电离层闪烁指数数据。2级数据产品包括：科学数据、图像产品和数据处理报告。

(1)科学数据

三频信标机2级科学数据文件包含以下内容：

◆VHF&UHF频段相对TEC值

◆L&UHF频段相对TEC值

◆VHF频段闪烁指数值

◆UHF频段闪烁指数值

◆L频段闪烁指数值

三频信标机2级数据格式见表5、表6所示。

表5三频信标机2级数据文件属性说明

表6三频信标机2级数据表格式说明

(2)2级图像产品

三频信标机2级图像产品包含以下5组曲线：

◆VHF&UHF频段相对TEC值-t

◆L&UHF频段相对TEC值-t

◆VHF频段闪烁指数-t

◆UHF频段闪烁指数-t

◆L频段闪烁指数-t

三频信标机2级图像产品示意图如图8所示。

(3)2级数据处理报告

三频信标机2级数据处理报告文件组成要素：

(1)处理软件版本号：V0.1

(2)起始时间：yyyymmdd HH:MM:SS.ZZZ

(3)输入数据：1级数据文件名

(4)辅助数据

(5)处理过程

①处理过程正常

②异常现象及可能的原因

(6)结束时间：yyyymmdd HH:MM:SS.ZZZ

(7)输出数据：2级数据文件名

2.2.4 3级数据产品

3级数据是在2级数据的基础上，进行反演处理生成一条观测链路的绝对TEC和二维剖面产品。根据电子密度二维分布结果，提取NmE、NmF1、NmF2异常。3级数据产品包括：科学数据及其图像和数据处理报告。

科学数据

3级数据产品格式详见表7、表8。

表7三频信标3级数据文件属性说明

表8三频信标3级数据表格式说明

(2)图像产品

3级数据产品示意图如图9、图10所示。

(2)3级数据处理报告

3级数据处理报告文件类型为TXT文件。组成要素：

(1)处理软件版本号：V0.1

(2)起始时间：yyyymmdd HH:MM:SS.ZZZ

(3)输入数据：(a)2级数据文件名；(b)前5条重访轨道的3级数据，列出其文件名

(4)辅助数据

①地震事件：该轨道时间之前1天内、100km范围6级以上地震记录

②磁情指数：该轨对应时间的Kp、Dst等

③IRI模型

(5)处理过程

①处理过程正常

②异常现象及可能的原因

(6)结束时间：yyyymmdd HH:MM:SS.ZZZ

(7)输出数据：3级数据文件名

3数据处理流程

地面三频信标接收机跟踪锁定星载三频信标信号后，输出三个频段信号的正交分量Q和同向分量I观测数据和信噪比数据(Level 0数据)。对三频信标接收机的I/Q观测数据分别进行差分相位计算和信号强度计算，可计算得到VHF&UHF、UHF&L的差分相位值和三个频段各自的信号强度值(Level 1数据)。对计算得到的差分相位值进行相位连接处理，每秒计算出一组相对电离层总电子含量TEC值；对计算得到的信号强度值进行幅度闪烁指数计算，每1秒计算出一组闪烁指数。进行时空匹配，计算得到各组数据在观测时间所对应的空间位置(Level 2数据)；通过对多个台站的观测数据进行集中处理，利用多站法每秒钟计算一个电离层绝对TEC数据，再利用层析成像算法对单个链路进程反演得到二维电子密度剖面数据，并提取NmF2异常(Level 3)。三频信标系统各级数据处理流程图见图11所示。

4数据处理方法

4.1 0级数据生成1级数据

三频信标0级数据生成1级数据的处理流程如图12所示。

4.1.1差分相位计算

第1步：读取0级VHF频段I路、Q路数据

三频信标0级数据主要是VHF与UHF、L与UHF频段两两组合，进行相位差分处理后得到的I路(同向)和Q路(正交)两组数据。其中UHF是接收机锁定信号时进行差分处理的基准频率，因此得到的有效差分数据只有VHF/UHF和L/UHF两种组合。每秒钟50组数据。

第2步：计算VHF/UHF差分相位

按照公式(1)进行反正切处理，计算得到VHF/UHF频段对应的差分相位。

Φ(t)＝arctan(Q/I) (1)

第3步：读取0级L频段I路、Q路数据

读取L频段I路和Q路数据，每秒钟50组数据。

第4步：计算L/UHF差分相位

按照公式(1)进行反正切计算，得到L/UHF频段对应的差分相位数据。

4.1.2信号强度计算

第1步：读取0级VHF、UHF、L频段的I路、Q路数据

分别读取三频信标0级VHF、UHF、L频段的I路、Q路数据。

第2步：计算VHF、UHF、L频段的信号强度

按照公式(2)进行信号强度计算。计算结果为功率值，分别得到VHF、UHF、L频段对应的三组信号强值数据，单位为dBm。

其中，10lg(I²+Q²)的处理是把接收机处理产出的相对电压值转换为功率值，231(dB)是通道增益量(包含天线增益、射频增益和算法增益)。

4.2 1级数据生成2级数据

三频信标1级数据生成2级数据的处理流程如图13所示。

4.2.1相对TEC计算

TEC是指电离层电子密度沿卫星-接收机链路的路径积分量。相干信标接收机所记录的差分多普勒相位限制在0～360°之间，因此必须将数据连接起来，才能变成相应的相对TEC。相位连接方法是在数据相位0～360°中取一个合适的门限值对相位进行判断和数据连接，具体实现步骤如图14所示。

具体处理方法流程为：

第1步：读取差分相位数据

分别读取VHF&UHF、UHF&L频段的差分相位数据。

第2步：取初始值，设置相位门限

取t＝1，k(1)＝1，D＝300。

第3步：差分计算

对于有n个点的相位数据，按照公式(3)所示对相位数据进行差分计算。

Δ(t)＝Φ(t+1)-Φ(t)，t＝1,2,....,n-1 (3)

第4步：相位翻转个数计算

将各数据点的Δ值与相位门限值D按照公式(4)的关系进行比较，得到各相位数据在连接时所取2π的翻转个数k。

第5步：相位连接处理

按照公式(5)依次对各数据点的差分相位序列进行相位值的增加和减少，即相位连接。

Φ'_c(t)＝Φ(t)+2kπ (5)

第6步：计算最小相位值

在连接完成后，按公式(6)计算相位序列中的最小相位值Z。

Z＝min(Φ'_c(t)) (6)

第7步：最小相位归零计算

Φ_c(t)＝Φ'_c(t)-Z (7)

其中Φ'_c(t)为连接起来的相位序列，按照公式(7)依次减掉最小相位值Z，最终得到连接完成的差分多普勒相位序列Φ_c(t)。

第8步：计算相关系数C_D

按照公式(8)计算相关系数C_D。其中，c＝299792458m/s，f_r＝16.668MHz，m_i表示倍频数，如果取VHF/UHF组合，则取m₁＝9，m₂＝24；如果取UHF/L组合，则取m₁＝24，m₂＝64。

第9步：计算相对TEC值

按照公式(4-9)可计算得到电离层相对TEC值。TEC值的单位为TECU(1TECU＝1016el.m-2)。

TEC_r＝C_D×Φ_c(t) (9)

4.2.2闪烁指数计算

三频信标机幅度闪烁指数S4的过程如图15所示。

第1步：读取信号强度值

分别读取VHF、UHF、L频段的信号强度值数据。

第2步：数据分组

以每50个数(1秒钟长度)为1组进行分组。

第3步：计算各组数据平方的均值

按照公式(10)计算得到每组数据各个数值的平方的平均值。

第4步：计算各组数据平均值的平方

按照公式(11)计算得到每组数据的平均值的平方。

第5步：计算闪烁指数值

按照公式(12)每1秒钟计算得到一个闪烁指数值。一般0.1≤S4<0.3被认为是电离层弱闪烁，0.3≤S4<0.6表示电离层中等闪烁，S4>0.6表示电离层强闪烁。

4.3 2级数据生成3级数据

三频信标2级数据生成3级数据的处理流程如图16所示。

(1)获取本站链对应的所有接收站的相对TEC数据，利用多站法，计算出每个接收站的积分常数，将积分常数加至接收机得到的的相对TEC，即得到各个接收站的绝对TEC。

绝对TEC具体计算步骤如下：

首先利用双频信标计算绝对TEC。

如图17，设接收机A到卫星S的斜TEC为I_s1，其投影的垂向TEC为I_v1，则二者具有如下关系：

其中，D₁称为路径S1上斜向全电子含量与垂直全电子含量的转换函数，χ_A为路径S1与接收站A垂线方向的夹角。卫星过境期间斜向全电子含量I_s1与差分多普勒相位关系式表示为：

其中I_A(t_k)为接收机A在t_k(k＝1，…，n)时刻观测给出的差分相对TEC；Φ₁(0)为待求接收站A的未知积分常数。其值等于1.3×10,3，其中c为光速，f₀＝16.668MHz，n₁＝3,n₂＝8；Φ_d1(t_k)是接收机A在各时刻观测对应的差分多普勒相位；公式(14)可表示为

其中，D_1k和I_v1k分别为第k个观测点与A对应射线的转换函数和垂直TEC。同理，对于接收站B有：

公式(15)、(16)可以变为

近似地认为，对于地面不同位置的两个接收站，两个接收站的垂向TEC相等。因此对于双站法I_v1k≈I_v2k，利用最小平方误差限定I_v1k和I_v2k两者间的误差，双站法的目标函数可表示下式：

通过对E求Φ₁(0)和Φ₂(0)的偏导数，再令其偏微分值为零，配合行列式求解。即联立方程 即可求得两站的未知积分常数Φ₁(0)和Φ₂(0)。

从而带入公式(14)可由观测得到的差分相对TEC得到各个站的双频绝对TEC如下式所示：

推广到m个接收站，也可按照相同的概念求解多个站的初始相位，即所谓的多站法(multi-station method)。多站法的目标函数为公式(4-20)。

通过求解E的极小值也可得到各个站的未知相位积分常数Φ_i(0)(i＝1,···,m)，从而求得各个站的双频绝对TEC。

②利用双频辅助三频求解三频信标绝对TEC

三频信标与双频信标相比多了一个载波两个频率的载波可做如下差分：

在三频信标中有:

f₁＝9f₀＝150.012MHz

f₂＝24f₀＝400.032MHz

f₃＝64f₀＝1066.752MHz

因此三组载波利用差分多普勒技术可以得到三个差分相位表达式(23)：

式中，f0＝16.668MHz，n1＝9,n2＝24,n3＝64，ΔP₁₂、ΔP₁₃、ΔP₂₃为差分相位的绝对值，包含接收机给出的差分多普勒相位和相位积分常数。

则：

由公式(24)可得：

ΔP₁₂(n₁ ²+n₂ ²)＝ΔP₁₃n₂ ² (25)

上式可改写为：

(K₁₂+Δφ₁₂)(n₁ ²+n₂ ²)＝(k₁₃+Δφ₁₃)n₂ ²或

K₁₂(n₁ ²+n₂ ²)-K₁₃n₂ ²＝Δφ₁₃n₂ ²-Δφ₁₂(n₁ ²+n₂ ²) (26)

其中，K₁₂和K₁₃是ΔP₁₂和ΔP₁₃整数部分，Δφ₁₂和Δφ₁₃是其小数部分，由接收机数据给出。基于素数原理，有：

x₁₂(n₁ ²+n₂ ²)-x₁₃n₂ ²＝1 (27)

其中，n₁和n₂互质，联立公式(26)、(27)可解出K₁₂和K₁₃，其表达式为：

K₁₂＝[Δφ₁₃n₂ ²-Δφ₁₂(n₁ ²+n₂ ²)]x₁₂+kn₂ ²

K₁₃＝[Δφ₁₃n₂ ²-Δφ₁₂(n₁ ²+n₂ ²)]x₁₃+k(n₁ ²+n₂ ²) (28)

其中k为一整数，则由K₁₂和K₁₃可求得ΔP₁₂和ΔP₁₃，再把ΔP₁₃或ΔP₁₂带入式(23)并整理可得到：

在三频信标中n₁＝3，n₂＝8，f₀＝16.668MHz。带入上式得到

TEC＝8.3165×10¹⁶×[(Δφ₁₃x₁₂-Δφ₁₂x₁₃)_mod1+k₂] (30)

其中8.3165×1016(单位：el/m2)称为三频信标测量TEC时的模糊系数，k₂为一正整数，它与双频信标的相位积分常数一样，为待定值。

由公式(19)、(30)，我们可以将双频和三频TEC的求解简化表示为：

TEC＝k₁×1.2995×10¹⁵+φ₁ (31)

TEC＝k₂×8.3165×10¹⁶+φ₂ (32)

其中，k₁和k₂分别表示双频和三频相位积分常数，其值均为整数，φ₁和φ₂分别表示双频和三频差分相位的小数部分。公式(31)通过多站法可以计算出包含误差的相位积分常数k₁，再结合地面接收机测得的差分相位小数部分φ₁，可以计算出包含误差的绝对TEC。由于同一反演区域同一时刻的TEC是相同的，将公式(31)计算得到的绝对TEC带入公式(32)，即得到k₂的计算公式：

因为公式(31)计算得到的绝对TEC存在误差，所以公式(33)计算得到的k₂一般也存在误差，我们需要用三频相位测量结果对其进行精确修正。研究表明，双频测量TEC的误差在1～3TECU量级(1TECU＝1.0×1016el/m2)，而三频相位模糊系数为8.3TECU(即以8.3TECU为单位分组)，结合k₂为整数的基本前提，对k₂进行精确求解的算法流程图如图18。

(2)将站链上空100km～500km范围按照高度间隔20km，纬度间隔0.5°进行间隔划分(图19)，在利用卫星信号观测进行的电离层层析成像过程中，所获得的电离层电子总含量TEC是电离层电子密度沿卫星至接收机射线路径上的积分，可表示为：

TEC＝∫_pN_eds (34)

其中Ne表示电子密度分布,p是接收机和卫星间的路径。

电子密度反演计算过程如下：

①先利用球谐分析和经验正交函数组合的函数基模型法进行电离层CT。球谐函数和经验正交函数选择用于电子密度的表示：

其中θ和分别代表纬度和经度,h代表高度.是归一化的勒让德函数,m,n,k是球谐函数和经验正交函数(EOF)的阶数，其中EOF通过经验模式或测量的电子密度数据通过奇异值分解得到。

根据卫星轨道和接收机位置确定反演区域后，将待反演区域内连续分布的电子密度离散化成N个不同网格，并假设每一个网格内电子密度为恒定不变的。然后确定反演区域内的有效射线，有效射线是指与反演划定区域上、下高度边界(与地球表面同心的球壳)的交点均落在反演区域内的射线。每一个网格内射线穿过的路径投影长度可按如下方法计算：首先将卫星轨道及地面接收站的球坐标转换成为直角坐标(X,Y,Z)，其转换公式为：

其中R＝R_e+h，R_e表示地球半径。然后利用直线方程将每一条射线通过各个网格上的交点在直角坐标系下的坐标求出。地面接收站与卫星间射线的直线方程可表示为：

其中(X_sat,Y_sat,Z_sat)与(X_sta,Y_sta,Z_sta)分别表示卫星与三频信标接收机(站)的位置坐标。(X,Y,Z)表示直线上的任何一点，K表示比例常数，对于每个固定的点(X,Y,Z)来说，K都是一个可以确定的常量。由上式我们可以得到：

若交点位于距地球中心点距离为R的面上，则R已知，因此有：

X²+Y²+Z²＝R² (39)

将上式代入其中，再由一元二次方程AK²+BK+C＝0求解得到K，代回上式中即可得到相交点的位置坐标，然后再计算

求得各网格间的投影长度Δl。所有射线在反演区域内各网格间的投影长度的集合即可表示为电离层电子密度反演的投影矩阵。确定投影矩阵以后，通过转换，电离层电子密度反演的问题可以简化为：

y＝Hx (41)

其中x是基函数的权重系数，绝对TEC可以用一列的y来表示.H是由经验正交函数和球谐函数及几何路径投影长度组成的算子。由于卫星信标多站法估计绝对TEC存在较大误差，因此我们对以上矩阵进行了修正，我们采用对应的一个接收机TEC相互差分的方法来多站法误差的影响，表示为:

b＝y-y₀＝(H-H₀)x＝Ax (42)

其中b是TEC数据的差分结果，y0是给定接收机的一个参考测量值.其中前向算子A包含了H和参考算子H0。采用截断奇异值分解(TSVD)正则化的方法求解矩阵方程，根据矩阵的奇异值分解理论，可将矩阵A分解为

式中：U＝(u₁,u₂,...,u_n)，u_i为左奇异值向量；V＝(v₁,v₂,...,v_n)，v_i为右奇异值向量；Σ＝(σ₁,σ₂,...,σ_n)，σ_i为奇异值，σ₁≥σ₂≥…≥σ_n≥0,k为截断系数。通过式(44)求解得到基函数权重x后，重新带入式(41)便可得到电子密度的初始分布值。

②将函数基模型的反演结果作为迭代初值，采用像素基模型的ART算法，对式(42)进行第二次迭代重构，以获取更小尺度的电离层扰动结构，计算方法如下：

式中，q表示迭代轮次，具体次序可按顺序或通过随机数的方法产生，i＝mod(q，L)+1,L为矩阵y的行数，t_i表示矩阵y的第i行元素，R_i表示R的第i个元素，λ为松弛因子。ART算法收敛较快，一般取迭代10～20轮左右即可，也可以在具体执行过程中设置适当的迭代终止阀值。通过以上两步的反演即可获得最终的电离层电子密度分布。

以上所述，仅是本发明的较佳实施例而已，并非对本发明作任何形式上的限制，本领域技术人员利用上述揭示的技术内容做出些许简单修改、等同变化或修饰，均落在本发明的保护范围内。