本文主要是介绍读取官方发布的DCB文件进行硬件延迟改正实现STEC提取，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！
一、引言
如果看到这里，我就将大家当成理解甚至已经实现STEC提取步骤，目前处于需要进行硬件延迟改正的阶段。但是如果大家还理不太清楚，可看参考我之前写过的STEC提取步骤，里面较为详细地推导了STEC提取方程，该篇博客其实专业性较强，需要大家之前阅读过一些相关文献或者是书籍，如果看不太懂话我写的或者想深入了解这部分理论的话，哈哈，不要紧！我个人建议大家去读武汉大学、中科院相关的博士论文或者高水平期刊，最好是外文的（绝大多中文期刊、学位论文里面的原理很多都含糊不清甚至有错误）。
硬件延迟改正也叫DCB改正，顾名思义就是改正伪码上的卫星硬件延迟和接收机硬件延迟，为什么不消除载波观测值的硬件延迟呢？从上篇博客STEC提取方程的推导中可以看到载波上的硬件延迟是被吸收进模糊度参数中一并估计的。所以在STEC提取中，最大的误差源就是伪码上的硬件延迟，若是不进行DCB改正，提取到的STEC会出现值异常。实验中我们可以观察到，未改正的STEC出现-10e+2TECU也属于正常现象。
由此可见，未进行DCB改正的STEC是无法用于电离层建模的，所以说正确估计卫星和测站的DCB是保证STEC提取精度的重要一环，下面介绍利用CAS发布的BSX文件进行DCB改正的方法。
二、DCB改正
2.1 DCB改正方式
DCB改正方式我大致归为两类：一是利用官方发布的DCB文件读取卫星和测站的DCB进行改正，二是自估计DCB，其中第一类是最为基础，直接的，但是精度不高（实验部分会进行对比），甚至有些测站不再文件测站列表中导致该测站上无法进行单站TEC提取，第二类是指采用各种算法估计测站和卫星的DCB，这里就包含单站估计以及多站全球建模两种方式。下面如何介绍利用CAS发布的BSX文件进行DCB改正。
2.2 读取BSX文件
首先我们需要了解BSX文件，我简单介绍一下吧怕有些同学不熟悉或者还没接触。
# 这里给出了提供改正信息的卫星系统，
# 以及该系统下卫星DCB的改正频率
- Overview of the DCBs included in the file GPS C1C-C1W, C2C-C2W, C2W-C2S, C2W-C2L, C2W-C2X GPS C1C-C2W, C1C_C5Q, C1C-C5X, C1W-C2W GLONASS C1C-C1P, C2C-C2P GLONASS C1C-C2C, C1C-C2P, C1P-C2P GALILEO C1C-C5Q, C1C-C7Q, C1C-C8Q, C1X-C5X, C1X-C7X, C1X-C8X BDS-2 C2I-C7I, C2I-C6I
-FILE/COMMENT然后是卫星硬件延迟，需要注意的是卫星硬件延迟是频率相关项，具体选择什么频率组合取决于你是用哪两个频率进行的STEC提取。
#
# GPS->G01->C1W-C2W组合上的卫星硬件延迟为-7.7370ns
# 由于卫星硬件延迟与频率相关，所以一定要根据你构建的STEC提取方程进行选择DSB G063 G01 C1W C2W 2017:244:00000 2017:245:00000 ns -7.7370 0.0290DSB G061 G02 C1W C2W 2017:244:00000 2017:245:00000 ns 9.2260 0.0290DSB G069 G03 C1W C2W 2017:244:00000 2017:245:00000 ns -5.4440 0.0345DSB G036 G04 C1W C2W 2017:244:00000 2017:245:00000 ns -0.5330 0.0300DSB G050 G05 C1W C2W 2017:244:00000 2017:245:00000 ns 2.8960 0.0300DSB G067 G06 C1W C2W 2017:244:00000 2017:245:00000 ns -7.0080 0.0290DSB G048 G07 C1W C2W 2017:244:00000 2017:245:00000 ns 3.3220 0.0360DSB G072 G08 C1W C2W 2017:244:00000 2017:245:00000 ns -7.5910 0.0290DSB G068 G09 C1W C2W 2017:244:00000 2017:245:00000 ns -4.8750 0.0290DSB G073 G10 C1W C2W 2017:244:00000 2017:245:00000 ns -5.4590 0.0290DSB G046 G11 C1W C2W 2017:244:00000 2017:245:00000 ns 3.8460 0.0290DSB G058 G12 C1W C2W 2017:244:00000 2017:245:00000 ns 3.8670 0.0290DSB G043 G13 C1W C2W 2017:244:00000 2017:245:00000 ns 3.1550 0.0290DSB G041 G14 C1W C2W 2017:244:00000 2017:245:00000 ns 2.0710 0.0295DSB G055 G15 C1W C2W 2017:244:00000 2017:245:00000 ns 3.0080 0.0290DSB G056 G16 C1W C2W 2017:244:00000 2017:245:00000 ns 2.9540 0.0290DSB G053 G17 C1W C2W 2017:244:00000 2017:245:00000 ns 3.3400 0.0290DSB G054 G18 C1W C2W 2017:244:00000 2017:245:00000 ns 3.1740 0.0290DSB G059 G19 C1W C2W 2017:244:00000 2017:245:00000 ns 6.0230 0.0290DSB G051 G20 C1W C2W 2017:244:00000 2017:245:00000 ns 1.4850 0.0290最后是测站DCB，该物理量同样需要选择相应的频率组合。
# ALBH->C1C-C1W
# 测站DCB为3.1370nsDSB G G ALBH C1C C1W 2017:244:00000 2017:245:00000 ns 3.1370 0.0165DSB G G ALG2 C1C C1W 2017:244:00000 2017:245:00000 ns 2.8910 0.0170DSB G G ALG3 C1C C1W 2017:244:00000 2017:245:00000 ns 2.7920 0.0170DSB G G ALGO C1C C1W 2017:244:00000 2017:245:00000 ns 1.4050 0.0170DSB G G AMC2 C1C C1W 2017:244:00000 2017:245:00000 ns 0.3450 0.0165DSB G G AREV C1C C1W 2017:244:00000 2017:245:00000 ns 1.4650 0.01652.3 读取代码以及示例
""" Read BSX file to get dcb of satellites and receivers """""" ----------------------------------------- """
import numpy as np""" ----------------------------------------- """""" --------------const---------------------- """
""" const_max_satellites_nums """
const_max_satellites_nums = 120*4*9
""" const_max_receiver_nums """
const_max_satellites_nums = 2400*9*4
""" ----------------------------------------- """""" class SatDCB """
class SatDCB:def __init__( self ):self.SVN = [ [0]*4 for i in range( const_max_satellites_nums ) ]self.PRN = [ [0]*3 for i in range( const_max_satellites_nums ) ]self.obs1 = [ [0]*3 for i in range( const_max_satellites_nums ) ]self.obs2 = [ [0]*3 for i in range( const_max_satellites_nums ) ]self.START = [ [0]*14 for i in range( const_max_satellites_nums ) ]self.END = [ [0]*14 for i in range( const_max_satellites_nums ) ]self.UNIT = [ [0]*4 for i in range( const_max_satellites_nums ) ]self.dcb = np.zeros( shape=(const_max_satellites_nums) )self.std = np.zeros( shape=(const_max_satellites_nums) )""" class RecDCB """
class RecDCB:def __init__( self ):self.STATION = [ [0]*4 for i in range( const_max_satellites_nums ) ]self.obs1 = [ [0]*3 for i in range( const_max_satellites_nums ) ]self.obs2 = [ [0]*3 for i in range( const_max_satellites_nums ) ]self.START = [ [0]*14 for i in range( const_max_satellites_nums ) ]self.END = [ [0]*14 for i in range( const_max_satellites_nums ) ]self.UNIT = [ [0]*4 for i in range( const_max_satellites_nums ) ]self.dcb = np.zeros( shape=(const_max_satellites_nums) )self.std = np.zeros( shape=(const_max_satellites_nums) )""" class ReadBSX """
class ReadBSX:def __init__( self, fname ):""" DCB filename """self.filename = fname""" DCB file object """self.dcbf = open( self.filename, "r" )""" satellite dcb """self.satDCB = SatDCB()""" receiver dcb """self.recDCB = RecDCB()""" fileline index(used to read bsx file) """self.sindex = 0self.rindex = 0def readdcb( self ):for line in self.dcbf:if "DSB"==line[1:4] and line[15:19].isspace():self.satDCB.SVN[self.sindex] = line[6:10]self.satDCB.PRN[self.sindex] = line[11:14]self.satDCB.obs1[self.sindex] = line[25:28]self.satDCB.obs2[self.sindex] = line[30:33]self.satDCB.START[self.sindex] = line[35:49]self.satDCB.END[self.sindex] = line[50:64]self.satDCB.UNIT[self.sindex] = line[65:69]self.satDCB.dcb[self.sindex] = line[69:91]self.satDCB.std[self.sindex] = line[91:103]self.sindex = self.sindex + 1elif "DSB"==line[1:4] and line[15:19]:self.recDCB.STATION[self.rindex] = line[15:19]self.recDCB.obs1[self.rindex] = line[25:28]self.recDCB.obs2[self.rindex] = line[30:33]self.recDCB.START[self.rindex] = line[35:49]self.recDCB.END[self.rindex] = line[50:64]self.recDCB.UNIT[self.rindex] = line[65:69]self.recDCB.dcb[self.rindex] = line[69:91]self.recDCB.std[self.rindex] = line[91:103]self.rindex = self.rindex + 1""" --------------interface----------------- """
""" get satellite dcb """
def get_sdcb( sdcb, prn, obst1, obst2 ):for i in range( const_max_satellites_nums ):if (sdcb.PRN[i]==prn) and (sdcb.obs1[i]==obst1) and (sdcb.obs2[i]==obst2):return sdcb.dcb[i]return -999
""" get receiver dcb """
def get_rdcb( rdcb, sta, obst1, obst2 ):for i in range( const_max_satellites_nums ):if rdcb.STATION[i]==sta and rdcb.obs1[i]==obst1 and rdcb.obs2[i]==obst2:return rdcb.dcb[i]return -999
""" ---------------------------------------- """
下面给出一个示例程序进行演示：
# 输入BSX文件路径（以MAXOS系统为例）
bsx = ReadBSX.ReadBSX("./201791/CAS0MGXRAP_20172440000_01D_01D_DCB.BSX")
# 进行卫星和测站DCB信息
bsx.readdcb()# 获取JFNG测站C1C-C2W上的DCB
rec_dcb = ReadBSX.get_rdcb(bsx.recDCB, “JFNG”, “C1C”, “C2W”)
# 存在BSX文件不提供该测站DCB的情况
if rec_dcb==-999:raise ValueError("该测站DCB无法从BSX文件中获取")# 获取G01卫星C1C-C2W上的DCB
sat_dcb = ReadBSX.get_sdcb(bsx.satDCB, “G01”, "C1C", "C2W")
# 出于程序完备性考虑
# 可能输入卫星prn或是输入该导航系统不存在的观测类型
# 或者你选择的双频组合不在文件提供的范围内
if sat_dcb==-999:raise ValueError("该卫星DCB无法从BSX文件中获取")三、利用CAS发布的BSX文件进行VTEC提取
3.1 VTEC提取结果对比
图1 利用CAS发布的DCB文件提取VTEC与自估计DCB的VTEC提取二者效果对比
3.2 总结
从图中不难看出利用官方发布的DCB进行TEC提取，在单站上的提取精度相对于自估计DCB的提取精度较差，有时甚至会发生从文件中获取不到对应测站DCB的情况，此时就必须采用算法进行DCB估计了。
到目前为止，已经有很多的DCB估计方式包括单站估计、多站建模。教大家一个最简单的实现方式也是麻省理工学院（MIT）2006年推出MAPGPS估计方式，卫星DCB我们暂且通过文件读取，测站DCB在MAPGPS估计中作为估计参数。
图2 我采用算法估计DCB后提取的GPS系统的VTEC
从上图可以看出，一般来说VTEC在夜间时表现平稳，在白天扰动较大，在图上表现为夜晚各卫星VTEC弧段最为贴合。而MIT在06年提出的MAPGPS方法正是基于这样的思想来寻找一个值，该值使得VTEC弧段变得最为聚合，该值即为对应测站DCB。
该方法有一定局限性，精度大多数情况不如其他常用估计方式（该方法的提取效果不如图2所展现的），但是实现上较为简单适合作为入门DCB估计算法，MIT2016年也就是10年后又推出了新的DCB估计方式以及统计框架，大家有兴趣可以去了解一下。（MIT也发布了电离层地图、视线向TEC等产品哟）
不知道我有没有讲清楚，如果大家有什么建议或者问题可以私信或者发在评论区，我不嫌麻烦😭，我有时间了等闲下来的时候肯定会及时回复。
这篇关于读取官方发布的DCB文件进行硬件延迟改正实现STEC提取的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！
