科创中国

一种大地测量非等距时序噪声的自适应差分估计方法

技术领域

本发明涉及一种大地测量技术，具体涉及一种非等距时序噪声的自适应差分估计方法。

背景技术

以全球卫星导航系统(GNSS)、多普勒无线电定位定位系统(DORIS)和卫星激光测距(SLR)、卫星测高(SA)等为代表的现代空间大地观测技术广泛应用于导航定位、卫星定轨、变形监测和参考框架维持等方面。在上述应用中通过观测值时序噪声估计确定随机模型能够有效提高应用精度。

常规时序噪声估计或随机模型估计大多集中在GNSS领域。如专利号为***的中国专利“一种北斗接收机非差观测值噪声评估方法”利用北斗三频组合观测值和差分算子进行噪声估计；专利号为***的专利申请“一种GPS时间序列噪声模型建立方法”、专利号为***的专利申请“区域 GPS基准站坐标时间序列的噪声模型获得方法”利用极大似然估计法进行GPS站坐标时序噪声估计。在上述研究中，由于观测值中含有一定的污染，必须进行剔除粗差，从而导致数据缺失，使观测值时序由等间隔时序变为非等间隔时序，此时必须对非等间隔时序进行数据插值，从而可能导致时间序列随机特性改变。

此外，常规差分方法存在如下不足：一、不能处理非等距时序，必须先通过插值方法进行填补，但插值算法可能会破坏原时间序列的随机特性；二、时间基函数阶数和差分算子长度为固定关系，差分算子不具有自适应性；三、差分算子的确定没有采用最优估计准则，不能保证最优性。

发明内容

发明目的：本发明的目的在于针对常规差分方法的不足，提供了一种简单可行、适用于全球导航卫星系统GNSS、多普勒定轨系统DORIS、卫星激光测距SLR、卫星测高 SA观测值和IGS站坐标五种类型的观测值时序，且能够处理不同移动窗口节点长度、非等距时序的大地测量非等距时序噪声的自适应差分估计方法。

技术方案：本发明提供了一种大地测量非等距时序噪声的自适应差分估计方法，包括以下步骤：

(1)读取不同空间大地测量技术观测值文件或IGS站坐标文件中的观测值数据，并生成观测值时间序列；

(2)选取时序移动窗口节点长度作为差分算子长度，确定观测值时序逼近的时间基函数；

(3)基于所选时间基函数构建差分算子估计方程，利用最小二乘方法和最小值搜索准则确定自适应差分算子，并通过差分算子矩阵计算差分序列；

(4)通过改变差分算子长度获得不同长度下的差分算子和相应的差分序列，利用峰度和偏态检验确定最佳的自适应差分算子及其差分序列；

(5)利用误差传播律进行观测值时序的噪声估计。

进一步，步骤(2)设移动窗口节点长度为K的第i个窗口内观测值序列的移动时间窗口为t∈[ti+1，ti+K]，选取时间基函数Fp(t)；

在移动时间窗口t∈[ti+1，ti+K]内，观测值时间序列y(t)可采用时间基函数Fp(t)表示为逼近函数形式：

式中，y(t)表示窗口内的观测值时序，Fj(t)表示j阶时间基函数，p表示时间基函数的最高阶数，βj表示Fj(t)的确定性系数，εt表示观测噪声；

进一步，所述时间基函数Fp(t)的确定包括以下步骤：

①对移动时间窗口t∈[ti+1，ti+K]标准化为[-1，1]，即

其中，窗口内节点序号k∈[1，K]，τi+k为ti+k的归一化时间；

②在以下两种时间基函数中任意选择一种作为时间基函数Fp(ti+k)，p满足 p≤K-2：

A)时间基函数Fp(ti+k)取p阶一般多项式Tp(τi+k)，其具体形式为：

B)时间基函数Fp(ti+k)取p阶Bernoulli多项式Bp(τi+k)，其具体形式为：

其中，组合数bq表示q阶Bemoulli数，其递推公式为：

进一步，步骤(3)设移动窗口节点长度为K的第i个窗口的差分算子 di＝[di，l …di，K]T，自适应差分算子的估计包括以下步骤：

①构建差分算子估计方程

令差分序列满足零均值期望特性，同时考虑非零基准向量ei＝[eil … eiK]，则可建立差分算子估计方程：

其中：

②求解差分算子

基于最小二乘准则求解自适应差分算子：

③差分算子估计方程的非零基准向量优化

基于最小值准则进行非零基准向量的优化，非零基准向量ei有两种生成方式，其一是直接采用随机向量，其二是对Ai的前p+1行矩阵进行SVD分解取V任一行向量；

选定其中一种生成方式产生K个以上的非零基准向量，利用(6)进行自适应差分算子求解并选取达最小值时的di为最优自适应差分算子；

④差分算子单位化

为避免每个窗口内的差分算子对观测噪声具有不同放大系数，对所求差分算子进行单位化，计算公式如下：

⑤构建差分序列的差分算子矩阵

按移动窗口法依次求解每个窗口的最优自适应差分算子并构成最优差分算子矩阵D，其构成方式如下

其中，Di表示由自适应差分算子和零元素组成的行向量，n为观测值时序长度，m为所得差分序列的长度，即m＝n-K+1；

⑥通过差分算子矩阵计算获得差分序列

根据式(8)自适应差分算子矩阵求得整个差分序列为：

其中，y为观测值时间序列；

因此，若采用逐元素计算方式，第i个窗口的差分值表示为：

其中，di，k表示第i个差分算子的第k项，k＝1，2，...，K。

进一步，步骤(4)逐渐增加差分算子长度即增大K的取值，得到不同长度下的差分序列，利用偏度和峰度进行检验：求取不同差分序列的偏度系数a1和峰度系数a2；选择偏度a1接近0、峰度a2接近3的差分算子长度作为最优差分算子长度，此时的最优差分算子为自适应差分算子，偏度和峰度系数的估计公式为：

其中，σ分别为差分序列的平均值和标准差。

进一步，步骤(5)顾及原始观测值时相关协因数阵Qy，则利用误差传播律估计观测值序列的噪声水平为：

式中，trace[·]表示矩阵迹运算；

若原始观测值相互独立或忽略原始观测值时相关系数，即Qy取单位阵，则上式简化为：

有益效果：本发明针对现有常规差分算子无法应用于测绘领域非等距时间序列噪声估计的问题，基于自适应移动时间窗口和最小二乘平差原理构建新的差分算子，有效解决了单个观测站及以上的非等距、非平稳观测时间序列噪声估计问题。所提供的自适应差分估计方法能够估计单个观测站、等距或非等距时序的噪声水平，首先克服了常规差分算子长度由差分阶数唯一确定的缺陷，使其适用于不同的移动窗口节点长度确保了新差分算子的多节点自适应性；其次通过最小二乘准则求解的自适应差分算子以及差分序列偏度和峰度检验保证了新差分算子的统计最优性。

附图说明

图1为本发明方法流程图；

图2为实施例1中2015年GPS某站观测值P1&P2测码误差估计结果；

图3为实施例2中2015年Jason-2卫星DORIS L1相位误差估计结果；

图4为实施例3中2014年Jason-2卫星SLR某站观测误差估计结果；

图5为实施例4中2013年Jason-2卫星SA观测误差估计结果；

图6为实施例5中2010-2015年某IGS站坐标三维误差估计结果。

具体实施方式

下面对本发明技术方案进行详细说明，但是本发明的保护范围不局限于所述实施例。

实施例：一种大地测量非等距时序噪声的自适应差分估计方法，如图1所示，按数据类型可分为GNSS观测值、DORIS观测值、SLR观测值、SA观测值、IGS站坐标五类。

实施例1：收集2015年年积日320的GNSS某站的观测值文件(GPS卫星PRN号1～31)，其时间采样率为10s，从中提取P1码和P2码观测值，按照下述方法分别获得 P1和P2测码误差估计结果。GNSS观测值时间序列的噪声估计步骤如下：

(1)从CDDIS(ftp：//***)数据中心下载观测值文件(O文件)，或者用户采集的GNSS观测值文件，读取观测值文件中的伪距观测值P和相位观测值L：

其中，ρ表示测站到卫星的几何距离，c表示光速，dt，dT分别表示接收机钟差和卫星钟差，dion、dtrop为电离层延迟和对流层延迟误差，λ为波长，N表示整周模糊度，εP，εL分别表示伪距和载波相位观测值的噪声；

根据仅有GNSS观测值信息、观测值+卫星精密轨道信息、双频以上观测值信息三种情况，对伪距P、载波相位L分别构造三种观测值时间序列：

其中，i，j表示不同的载波频率；

(2)利用式(2)对移动时间窗口[ti+1，ti+K]标准化为[-1，1]，在一般多项式(3)、Bernoulli 多项式(4)中任选一种时间基函数；

(3)确定时间基函数最高阶数p，按照上述给出伪距P、载波相位L的各三种观测时间序列，阶数分别采用p＝4，p＝2，p＝2；

(4)选取差分算子长度K(移动窗口节点长度)的初值K＝p+2；

(5)根据式(5)构建自适应差分算子，利用式(6)求解差分算子；

(6)根据最小值准则搜索最优的自适应差分算子，并利用式(7)对所求自适应差分算子进行单位化；

(7)利用移动窗口法，依次求得各个单位化自适应差分算子，并根据式(8)构建自适应差分算子矩阵；

(8)利用式(9)求得观测值时序的差分序列；

(9)利用式(11)进行偏度α1和峰度α2检验；

(10)增加差分算子长度(移动窗口长度)K，求取不同长度下差分序列l的偏度和峰度系数；

(11)选择偏度接近0、峰度接近3的差分算子矩阵作为最终差分算子矩阵，该差分算子计算所得的差分序列为最终的差分序列l；

(12)利用式(12)～(13)估计GNSS伪距或相位观测值时序的噪声水平。

从图2可以看出，P1和P2码误差随历元变化平稳，没有出现趋势性变化，且误差水平与接收机噪声标称水平一致(分别约***、***)，表明本发明的方法正确有效。

实施例2：收集2015年3月1日至30日的Jason-2卫星DORIS观测值文件，其时间采样间隔为非等距(3s和7s交替采样)，从中提取L1载波观测值，按照下述方法获得L1相位误差估计结果。DORIS观测值时间序列的噪声估计步骤如下：

(1)从IDS(ftp：//***)数据中心下载观测值文件，或者用户收集的DORIS观测值文件，读取观测值文件中的载波相位观测值LD：

LD＝R+c(dt-dT)+λN-dion+dtrop+εD

其中，R表示DORIS地面站到卫星的几何距离，c表示光速，dt，dT表示接收机钟差和卫星钟差，dion为电离层延迟误差，dtrop为对流层延迟误差，λ为波长，N为整周模糊度，εD分别表示载波相位观测噪声；

根据仅有DORIS观测值信息、观测值+卫星精密轨道信息、双频以上观测值信息三种情况，对DORIS载波相位构造三种观测值时间序列：

其中，LD1，LD2为不同频率的载波相位观测值；

(2)利用式(2)对移动时间窗口[ti+1，ti+K]标准化为[-1，1]，在一般多项式(3)、Bernoulli 多项式(4)中任选一种时间基函数；

(3)确定时间基函数最高阶数p，按照上述给出载波相位LD的三种观测时间序列，阶数分别采用p＝4，p＝2，p＝2；

(4)选取差分算子长度K(移动窗口节点长度)的初值K＝p+2；

(5)根据式(5)构建自适应差分算子，利用式(6)求解差分算子；

(6)根据最小值准则搜索最优的自适应差分算子，并利用式(7)对所求自适应差分算子进行单位化；

(7)利用移动窗口法，依次求得各个单位化自适应差分算子，并根据式(8)构建自适应差分算子矩阵；

(8)利用式(9)求得观测值时序的差分序列；

(9)利用式(11)进行偏度α1和峰度α2检验；

(10)增加差分算子长度(移动窗口节点长度)K，求取不同长度下差分序列l的偏度和峰度系数；

(11)选择偏度接近0、峰度接近3的差分算子矩阵作为最终的差分算子矩阵，该差分算子计算所得的差分序列为最终的差分序列l；

(12)利用式(12)～(13)估计DORIS相位观测值时序的噪声水平。

从图3可以看出，L1相位误差随历元变化平稳，没有出现趋势性变化，且误差水平与DORIS站观测噪声评估水平一致(约***周)，表明本发明的方法正确有效。

实施例3：收集2014年1月1日至4月30目的Jason-2卫星SLR观测值文件(CRD 格式)，其时间采样间隔为非等距(数秒、数十秒不等)，从中提取激光测距观测值，按照下述方法获得激光测距误差估计结果。SLR观测值时间序列的噪声估计步骤如下：

(1)从EDC(ftp：//***)数据中心下载观测值文件，或者用户收集的SLR观测值文件，读取观测值文件中的距离观测值S：

其中，c表示光速，Δt表示信号传播时间，Δs为对流层、相对论、卫星质心和测站系统误差等，εS分别表示观测值噪声；

根据仅有SLR观测值信息、观测值+卫星精密轨道信息两种情况，对SLR激光观测值构造两种观测值时间序列：

(2)利用式(2)对移动时间窗口[ti+1，ti+K]标准化为[-1，1]，在一般多项式(3)、Bernoulli 多项式(4)中任选一种时间基函数；

(3)确定时间基函数最高阶数p，按照上述给出激光观测S的两种观测时间序列，阶数分别采用p＝4，p＝2；

(4)选取差分算子长度K(移动窗口节点长度)的初值K＝p+2；

(5)根据式(5)构建自适应差分算子，利用式(6)求解差分算子；

(6)根据最小值准则搜索最优的自适应差分算子，并利用式(7)对所求自适应差分算子进行单位化；

(7)利用移动窗口法，依次求得各个单位化自适应差分算子，并根据式(8)构建自适应差分算子矩阵；

(8)利用式(9)求得观测值时序的差分序列；

(9)利用式(11)进行偏度α1和峰度α2检验；

(10)增加差分算子长度K(移动窗口节点长度)，求取不同长度下差分序列l的偏度和峰度系数；

(11)选择偏度接近0、峰度接近3的差分算子矩阵作为最终的差分算子矩阵，该差分算子计算所得的差分序列为最终的差分序列l；

(12)利用式(12)～(13)估计SLR激光观测值时序的噪声水平。

从图4可以看出，SLR误差随历元变化平稳，没有出现趋势性变化，且误差水平与SLR测站噪声评估水平一致(约***)，表明本发明的方法正确有效。

实施例4：收集2013年年积日1-年积日254的Jason-2卫星SA观测值文件，其时间采样间隔为非等距(标称采样率1s，但并不严格)，从中提取测高观测值，按照下述方法获得SA测高误差估计结果。SA观测值时间序列的噪声估计步骤如下：

(1)从法国国家空间研究中心卫星海洋学存档数据中心(AVISO， ftp：//***)下载Jason-3数据，或者用户收集的其他SA观测值文件，读取观测值文件中测高观测值ha：

式中，h表示卫星的地心向距，h0表示卫星星下点地心向距，hN分别表示海水面到参考椭球面的距离，表示地心向距h和星下点的地心向距h0的不共线改正，εa分别表示观测噪声；

根据仅有SA观测值信息、观测值+卫星精密轨道信息、双频以上观测值信息三种情况，对SA观测值构造三种观测值时间序列：

其中，ha1，ha2为双频卫星测高观测值；

(2)利用式(2)对移动时间窗口[ti+1，ti+K]标准化为[-1，1]，在一般多项式(3)、Bernoulli 多项式(4)中任选一种时间基函数；

(3)确定时间基函数最高阶数p，按照上述给出测高观测值ha的三种观测时间序列，阶数分别采用p＝4，p＝2，p＝2；

(4)选取差分算子长度K(移动窗口节点长度)的初值K＝p+2；

(5)根据式(5)构建自适应差分算子，利用式(6)求解差分算子；

(6)根据最小值准则搜索最优的自适应差分算子，并利用式(7)对所求自适应差分算子进行单位化；

(7)利用移动窗口法，依次求得各个单位化自适应差分算子，并根据式(8)构建自适应差分算子矩阵；

(8)利用式(9)求得观测值时序的差分序列；

(9)利用式(11)进行偏度α1和峰度α2检验；

(10)增加差分算子长度(移动窗口节点长度)K，求取不同长度下差分序列l的偏度和峰度系数；

(11)选择偏度接近0、峰度接近3的差分算子矩阵作为最终的差分算子矩阵，该差分算子计算所得的差分序列为最终的差分序列l；

(12)利用式(12)～(13)估计卫星测高SA观测值时序的噪声水平。

从图5可以看出，SA误差随历元变化平稳，没有出现趋势性变化，且误差水平与测高标称水平一致(约6cm)，表明本发明的方法正确有效。

实施例5：收集了2010-2015年共5年的某IGS站三维坐标时间序列数据文件，其时间采样间隔为等间隔(天)，经过粗差剔除和阶跃项修复等预处理后，按照下述方法获得三维坐标误差估计结果。IGS站三维坐标时序的噪声估计步骤如下：

(1)从Sopac(http：//***)数据中心下载站点坐标文件，或者用户解算的站坐标时序文件，读取坐标观测值文件中三维方向坐标序列(X，Y，Z)：

式中，a，b分别表示初始位置和线性速度，N为周期项个数，Tj为周期，sj，cj为周期项系数，ε分别表示观测值噪声；

利用IGS站三维坐标序列构建观测值序列：

(2)利用式(2)对移动时间窗口[ti+1，ti+K]标准化为[-1，1]，在一般多项式(3)、Bernoulli 多项式(4)中任选一种时间基函数；

(3)确定时间基函数最高阶数p，按照上述给出IGS站坐标的三维时间序列，阶数分别采用p＝2，p＝2，p＝2；

(4)选取差分算子长度K(移动窗口节点长度)的初值K＝p+2；

(5)根据式(5)构建自适应差分算子，利用式(6)求解差分算子；

(6)根据最小值准则搜索最优的自适应差分算子，并利用式(7)对所求自适应差分算子进行单位化；

(7)利用移动窗口法，依次求得各个单位化自适应差分算子，并根据式(8)构建自适应差分算子矩阵；

(8)利用式(9)求得观测值时序的差分序列；

(9)利用式(11)进行偏度α1和峰度α2检验；

(10)增加差分算子长度(移动窗口节点长度)K，求取不同长度下差分序列l的偏度和峰度系数；

(11)选择偏度接近0、峰度接近3的差分算子矩阵作为最终的差分算子矩阵，该差分算子计算所得的差分序列为最终的差分序列l；

(12)利用式(12)～(13)估计观测值时序的噪声水平。

从图6可以看出，三维坐标误差随历元变化平稳，没有出现趋势性变化，且误差水平与测站标称水平一致(分别约为***，***，***)，表明本发明的方法正确有效。