一种区域高分辨率电离层TEC地图重构方法与流程

一种区域高分辨率电离层tec地图重构方法
技术领域
1.本发明属于空间天气现报预报领域，特别涉及该领域中的一种区域高分辨率电离层tec地图重构方法。


背景技术：

2.电离层是卫星导航、通信、雷达等诸多无线电信息系统最主要的误差来源之一。总电子含量(total electron content,tec)是表征电离层等离子体数密度的一个重要参数，掌握其时空变化特征和规律，对提升无线电信息系统性、提升空间科学研究及空间天气预报能力具有重要的意义。
3.随着卫星技术的发展，地基gnss已经成为全球tec地图构建最为重要的数量来源。自1998年以来，国际gnss服务(igs)开始公开发布全球tec地图(global ionospheric map,gim)数据。随着gnss站网数量的增加及数据处理技术的发展，目前已有多种生成gim的建模策略，包括球谐函数法、层析成像(ct)法、kalman滤波同化等。
4.随着全球范围内密集的gnss台站的建立，获取高分辨率的区域tec地图成为可能。近年来，美国mit利用全球超过2000个gnss监测站，发布了1
°×1°×
5分钟分辨率的散点tec数据产品；比利时皇家天文台(rob)利用欧洲区域euref密集的gnss监测站，实现了0.5
°×
0.5
°×
15分钟欧洲区域近实时tec地图产品的在线发布(ftp://gnss.oma.be/gnss/products/ionex/)。


技术实现要素：

5.本发明所要解决的技术问题就是针对区域精细化电离层监测和中小尺度扰动分析的需求，提供一种区域高分辨率电离层tec地图重构方法。
6.本发明采用如下技术方案：
7.一种区域高分辨率电离层tec地图重构方法，其改进之处在于，包括如下步骤：
8.步骤a，地基gnss观测数据下载与预处理：
9.步骤a1，下载gnss卫星观测数据和卫星星历数据文件；
10.步骤a2，读取数据文件，获取gnss原始双频观测数据，包含载波相位观测数据、伪距观测数据以及导航卫星星历；
11.步骤a3，对文件中gps的观测数据进行周跳和异常值的探测，并采用双频的载波相位观测数据平滑码伪距，最后输出平滑处理后的观测数据文件；
12.步骤b，地基gnss倾斜tec解算：
13.步骤b1，基于双频载波相位观测量，计算载波相位电离层延迟观测量具体方法如下：
14.[0015][0016][0017][0018][0019]
其中，i和j分别表示接收机和卫星的编号；c表示真空中的光速；和分别是卫星和接收机在频率f1和f2上的硬件延迟；λ1和λ2分别表示频率f1和f2对应的波长；和分别表示频率f1和f2上载波相位的模糊度；bi和bj分别是接收机和卫星的频间偏差，α为恒定常数；
[0020]
步骤b2，基于双频伪距观测量，计算伪距电离层延迟观测量具体方法如下：
[0021][0022]
步骤b3，将一连续弧段内某卫星的原始观测数据分为n组，根据步骤b1和b2分别获得n组载波相位和伪距电离层延迟观测量，计算伪距和载波相位电离层延迟观测量之和的平均值，具体表示为：
[0023][0024]
步骤b4，将步骤b3中获得的载波相位和伪距电离层延迟观测量之和的平均值回带至载波相位电离层延迟观测量中，获得卫星视线方向电离层延迟观测量具体表示为：
[0025][0026]
步骤b5，重复步骤b1至b4，逐卫星求解各基准站卫星视线方向上的高精度的电离层延迟观测量；
[0027]
步骤b6，采用ciraolo l提出的方法，分别计算出bi和bj，进而计算电离层倾斜tec，计算方法如下：
[0028][0029]
步骤c，区域多台站gnss垂直tec投影：
[0030]
步骤c1，将同一时段内的倾斜tec按照对应的穿刺点位置投影为垂直tec值：
[0031]
vtec＝tec/m(e,h)
[0032][0033]
其中，vtec表示接收机与卫星链路上投影得到的垂直tec，m(e,h)表示倾斜tec与垂直tec之间的投影函数，e为仰角，为接收机i和卫星j之间的仰角，re表示地球半径，h表示电离层薄层高度；
[0034]
步骤c2，按经度为区间进行划分，经度间隔取1
°
，将落在同一区域内的垂直tec数据进行分组，每组tec表示为iv；
[0035]
步骤c3，对相同经度区间内同组的垂直tec随纬度的变化进行线性拟合，采用阈值判定法对曲线中相比其它数据点偏差过大的点对应的tec进行滤除，判决条件如下：
[0036]
|iv
k-μ|＞2σ,k＝1,2,...,k
[0037]
其中，k表示数据编号，k表示iv的样本数目，μ为数据序列均值，σ为数据标准差，如满足该判决条件，则认为该tec值存在粗差，该数据被剔除；
[0038]
步骤c4，重复步骤c3，直至所有数据点都在可接受的范围内，以保证所有分组后倾斜tec转换得到垂直tec拟合的曲线都平滑为止；
[0039]
步骤d，区域高分辨率电离层tec地图重构：
[0040]
步骤d1，设定x代表存在空白间隙的vtec时空数据集，w是与x维度大小相同的数组，用于标识x值是否存在缺失；
[0041]
步骤d2，采用离散余弦变换—惩罚最小二乘回归算法寻求最小化的代价函数f表示如下：
[0042][0043]
其中|| ||是欧几里得范数，和分别代表拉普拉斯算子和schur乘积，s是控制平滑程度的标量；
[0044]
步骤d3，区域tec地图通过步骤d2及逆离散余弦变换来实现，表示为：
[0045][0046]
其中，γ是一个滤波张量，具体表示为：
[0047][0048]
其中，p
t
表示沿第t维的第p个元素，n
t
表示沿该维度x的大小。
[0049]
本发明的有益效果是：
[0050]
本发明所公开的方法，利用区域较为密集的地基gnss监测数据，基于离散余弦变换—惩罚最小二乘回归(dct—pls)算法，可以实现区域非常高的时间和空间分辨率tec地图重构。本发明方法可为实现区域高分辨率电离层监测及扰动分析提供技术支撑。
附图说明
[0051]
图1是本发明方法的流程示意图；
[0052]
图2是区域多台站gnss垂直tec投影效果图；
[0053]
图3是区域高分辨率电离层tec地图重构效果图。
具体实施方式
[0054]
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图和实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
[0055]
地震、海啸、台风、火箭发射等事件可能会导致小而明显的电离层扰动；在太阳耀斑或日冕物质抛射(cme)等空间天气事件中，短时间内电离层剧烈的扰动变化也时常发生。更小尺度电离层扰动的分析和建模亟需高时间和空间分辨率的电离层tec地图产品。目前大多数gim数据产品时间分辨率通常为1—2小时，地理纬度和经度网格大小为2.5
°×5°
，这种分辨率的tec地图无法满足精细化电离层监测及中小尺度扰动分析的需要。
[0056]
实施例1，本实施例公开了一种区域高分辨率电离层tec地图重构方法，如图1所示，包括如下步骤：
[0057]
步骤a，地基gnss观测数据下载与预处理：
[0058]
步骤a1，下载gnss卫星观测数据和卫星星历数据文件；
[0059]
步骤a2，读取数据文件，获取gnss原始双频观测数据，包含载波相位观测数据、伪距观测数据以及导航卫星星历；
[0060]
步骤a3，采用bernese软件对rinex格式文件中gps的观测数据进行周跳和异常值的探测，并采用双频的载波相位观测数据平滑码伪距，最后平滑处理后的观测数据仍以rinex格式的文件输出；
[0061]
步骤b，地基gnss倾斜tec解算：
[0062]
步骤b1，基于双频载波相位观测量，计算载波相位电离层延迟观测量具体方法如下：
[0063][0064][0065][0066][0067][0068]
其中，i和j分别表示接收机和卫星的编号；c表示真空中的光速；和分别是卫星和接收机在频率f1和f2上的硬件延迟；λ1和λ2分别表示频率f1和f2对应的波长；和分别表示频率f1和f2上载波相位的模糊度；bi和bj分别是接收机和卫星的频间偏差，α为恒定常数，一般可取值40.31；
[0069]
步骤b2，基于双频伪距观测量，计算伪距电离层延迟观测量具体方法如下：
[0070][0071]
其中，a，bi和bj的含义与步骤b1中相同；
[0072]
步骤b3，将一连续弧段内某卫星的原始观测数据分为n组，根据步骤b1和b2分别获得n组载波相位和伪距电离层延迟观测量，计算伪距和载波相位电离层延迟观测量之和的平均值，具体表示为：
[0073][0074]
上式中各变量的含义与步骤b1、b2中相同；
[0075]
步骤b4，将步骤b3中获得的载波相位和伪距电离层延迟观测量之和的平均值回带至载波相位电离层延迟观测量中，即可获得卫星视线方向电离层延迟观测量具体表示为：
[0076][0077]
步骤b5，重复步骤b1至b4，逐卫星求解各基准站卫星视线方向上的高精度的电离层延迟观测量；
[0078]
步骤b6，采用ciraolo l提出的方法(具体参考ciraolo l,azpilicueta f,brunini c,et al.calibration errors on experimental slant total electron content(tec)determined with gps[j].j.geod.,2007,81(2):111
–
120)，分别计算出bi和bj，进而计算电离层倾斜tec，计算方法如下：
[0079][0080]
步骤c，区域多台站gnss垂直tec投影：
[0081]
步骤c1，将同一时段内的倾斜tec按照对应的穿刺点位置投影为垂直tec值：
[0082]
vtec＝tec/m(e,h)
[0083][0084]
其中，vtec表示接收机与卫星链路上投影得到的垂直tec，m(e,h)表示倾斜tec与垂直tec之间的投影函数，e为仰角，为接收机i和卫星j之间的仰角，re表示地球半径，h表示电离层薄层高度，一般取350千米；
[0085]
步骤c2，按经度为区间进行划分，经度间隔取1
°
，将落在同一区域内的垂直tec数据进行分组，每组tec表示为iv；
[0086]
步骤c3，对相同经度区间内同组的垂直tec随纬度的变化进行线性拟合，采用阈值判定法对曲线中相比其它数据点偏差过大的点对应的tec进行滤除，判决条件如下：
[0087]
|iv
k-μ|＞2σ,k＝1,2,...,k
[0088]
其中，k表示数据编号，k表示iv的样本数目，μ为数据序列均值，σ为数据标准差，，如满足该判决条件，则认为该tec值存在粗差，该数据被剔除；
[0089]
步骤c4，重复步骤c3，直至所有数据点都在可接受的范围内，以保证所有分组后倾斜tec转换得到垂直tec拟合的曲线都较为平滑为止，处理结果如图2所示。
[0090]
步骤d，区域高分辨率电离层tec地图重构：
[0091]
步骤d1，设定x代表存在空白间隙的vtec时空数据集，w是与x维度大小相同的数组，用于标识x值是否存在缺失；
[0092]
步骤d2，采用离散余弦变换—惩罚最小二乘回归算法(dct—pls)寻求最小化的
实现tec地图的重构，代价函数f表示如下：
[0093][0094]
其中|| ||是欧几里得范数，和分别代表拉普拉斯(laplace)算子和schur乘积，s是控制平滑程度的标量；
[0095]
步骤d3，区域tec地图通过步骤d2及逆离散余弦变换(idct)来实现，表示为：
[0096][0097]
其中，γ是一个滤波张量，具体表示为：
[0098][0099]
其中，p
t
表示沿第t维的第p个元素，n
t
表示沿该维度x的大小。图3所示为区域高分辨率电离层tec地图重构效果图。