CN103148813B - 用于处理gps变形监测数据的方法 - Google Patents

Abstract

本发明提出了一种用于处理GPS变形监测数据的方法，包括：第一步：获取基准站和监测站当前历元观测数据及广播星历数据，利用相位观测值组成双差观测方程；第二步：设置两历元间的最大变形量；第三步：采用Cholesky分解构建模糊度搜索空间；第四步：搜索模糊度，并且使方差比大于3；第五步：获取监测点的单历元整周模糊度固定解坐标X、Y、Z和其协方差矩阵；第六步：利用平均间隙法对当前历元和前一历元间的变形值进行变形检验；第七步：经第六步的检验后，将两个历元观测方程叠加的抗差最小二乘静态解作为当前历元的结果或者根据抗差序贯平差获取当前历元的结果。本方法很准确。

Description

用于处理GPS变形监测数据的方法

技术领域

本发明涉及一种用于处理数据的方法，特别涉及一种用于处理GPS变形监测数据的方法。

背景技术

随着GPS技术的发展，GPS已越来越多地用于变形监测。但在变形监测环境中，通常多路径或衍射信号较强，导致周跳发生频率较高，而常用的最小二乘搜索法（LS）、快速模糊度搜索法（FARA）、降低模糊度之间相关性的最小二乘非相关平差法（LAMBDA）、Cholesky分解搜索法等都需要利用一段时间内的多个历元的观测数据信息，且这多个历元的观测数据不能出现周跳，否则将很难确定整周模糊度。

因此在变形监测应用中，GPS数据的处理一般采用单历元解算方法以避开周跳探测，即根据监测点位置变化不大的特点，利用变形量约束进行单历元模糊度求解。同样由于监测环境的限制，可视GPS卫星通常较少，导致卫星分布图形结构强度较弱，加上多路径效应的影响，使得所获取的结果精度不高，在变形序列中会出现跳变的情况，从而难以根据监测结果进行预警。

针对以上现有技术中的用于处理GPS变形监测数据的方法的不足，研究精度和可靠性都得到提升的新的GPS数据处理方法，对变形监测技术的发展和其在工程中的实际应用具有积极的意义。

发明内容

针对现有技术中的用于处理GPS变形监测数据的方法无法回避跳变，导致难以根据监测结果精确预警的不足，本发明提出了一种新颖的用于处理GPS变形监测数据的方法。

本发明提出了一种用于处理GPS变形监测数据的方法，包括：第一步：获取基准站和监测站当前历元观测数据及广播星历数据，利用相位观测值组成双差观测方程；第二步：设置两历元间的最大变形量；第三步：采用Cholesky分解构建模糊度搜索空间；第四步：搜索模糊度，并且使方差比大于3；第五步：获取监测点的单历元整周模糊度固定解坐标X、Y、Z和其协方差矩阵；第六步：利用平均间隙法对当前历元和前一历元间的变形值进行变形检验；第七步：经第六步的检验后，如果两历元间未发生位移，则将两个历元观测方程叠加的抗差最小二乘静态解作为当前历元的结果；若第六步的检验结果为两历元间有位移发生，则根据抗差序贯平差获取当前历元的结果。

在一个实施例中，在第三步中：

$-\frac{{\mathrm{\δd}}_{\max }}{{\mathrm{\λl}}_{11}}\≤{\mathrm{\δ}}_{\▿\mathrm{\Δ}{N}_1}\≤\frac{{\mathrm{\δd}}_{\max }}{{\mathrm{\λl}}_{11}}$

$\frac{-\sqrt{{{\mathrm{\δd}}_{\max }}^2-B^2}-l_{21}\mathrm{\λ}{\mathrm{\δ}}_{\▿\mathrm{\Δ}{N}_1}}{{\mathrm{\λl}}_{22}}\≤{\mathrm{\δ}}_{\▿\mathrm{\Δ}{N}_2}\≤\frac{\sqrt{{{\mathrm{\δd}}_{\max }}^2-B^2}-l_{21}\mathrm{\λ}{\mathrm{\δ}}_{\▿\mathrm{\Δ}{N}_1}}{{\mathrm{\λl}}_{22}}$

$\frac{-\sqrt{{{\mathrm{\δd}}_{\max }}^2-B^2-C^2}-l_{31}\mathrm{\λ}{\mathrm{\δ}}_{\▿\mathrm{\Δ}{N}_1}-l_{32}\mathrm{\λ}{\mathrm{\δ}}_{\▿\mathrm{\Δ}{N}_2}}{\mathrm{\λ}{l}_{33}}\≤{\mathrm{\δ}}_{\▿\mathrm{\Δ}{N}_3}\≤\frac{\sqrt{{{\mathrm{\δd}}_{\max }}^2-B^2-C^2}-l_{31}\mathrm{\λ}{\mathrm{\δ}}_{\▿\mathrm{\Δ}{N}_1}-l_{32}\mathrm{\λ}{\mathrm{\δ}}_{\▿\mathrm{\Δ}{N}_2}}{{\mathrm{\λl}}_{33}}$

其中：

$B=l_{11}\·\mathrm{\λ}{\mathrm{\δ}}_{\▿\mathrm{\Δ}{N}_1,}$ $C=l_{21}\·\mathrm{\λ}{\mathrm{\δ}}_{\▿\mathrm{\Δ}{N}_1}+l_{22}\·\mathrm{\λ}{\mathrm{\δ}}_{\▿\mathrm{\Δ}{N}_2},$

δd_max为监测点的最大变形量，为双差整周模糊度中误差，λ为

载波相位的波长，

$L^{-1}=\left[\begin{array}{ccc}{l}_{11}& & \\ l_{21}& l_{22}& \\ l_{31}& l_{32}& l_{33}\end{array}\right]$ 为求解条件方程组和误差方程组所组成的对称正定的方程组矩阵AA^T的Cholesky分解的下三角矩阵。

在一个实施例中，在第六步中还包括如下子步骤：

第一子步骤：通过两历元观测值改正数获取经验方差：

$S_0^2=\frac{{\left(V^T\mathrm{PV}\right)}_I+{\left(V^T\mathrm{PV}\right)}_{\mathrm{II}}}{f}$

其中f为两历元自由度之和，即f＝n₁-3+n₂-3，其中n₁,n₂为两历元双差观测方程的个数，V为观测值残差向量，P为观测值权阵；

第二子步骤：由两历元的坐标差d_i(i＝1，2，3)获取单位权经验方差

${\stackrel{\‾}{S}}_0^2=\frac{d^T{P}_{d}d}{3}$

其中，d＝X_II-X_I，表示监测点在前后两历元的坐标之差向量：

$d=\left(\begin{array}{cc}-1& 1\end{array}\right)\left(\begin{array}{c}{X}_I\\ X_{\mathrm{II}}\end{array}\right)$

其中P_d为d的权阵：

$P_d=Q_d^{-1}$

此处Q_d由 $d=\left(\begin{array}{cc}-1& 1\end{array}\right)\left(\begin{array}{c}{X}_I\\ X_{\mathrm{II}}\end{array}\right)$ 根据协方差传播定律得到：

$Q_d=\left(\begin{array}{cc}-1& 1\end{array}\right)\left[\begin{array}{cc}{Q}_{X_I{X}_I}& 0\\ 0& Q_{X_{\mathrm{II}}{X}_{\mathrm{II}}}\end{array}\right]\left(\begin{array}{c}-1\\ 1\end{array}\right)=Q_{X_I{X}_I}+Q_{X_{\mathrm{II}}{X}_{\mathrm{II}}}$

第三子步骤：构建统计量：

$F_{h,f}=\frac{{\stackrel{\‾}{S}}_0^2}{S_0^2}=\frac{d^T{P}_{d}d}{h{S}_0^2}$

第四子步骤：选用显著水平α=0.01，将获得的F_h,f与从F分布表中查出的F_α(h,f)分位值比较，如果

F_h,f＞F_α(h,f)

则表明大于即有：

P{[F_h,f＞F_α(h,f)]}＝α

从而判断出两历元间有位移发生；

如果无法得出F_h,f＞F_α(h,f)，则判断出两历元间没有发生位移。

在一个实施例中，在第四步中，利用最小二乘法来搜索模糊度。

在一个实施例中，在第五步中，利用最小二乘法来获取监测点的单历元整周模糊度固定解坐标X、Y、Z和其协方差矩阵。

根据本发明的用于处理GPS变形监测数据的方法相对现有技术带来了如下进步：避免了周跳探测；采用Cholesky分解构建模糊度搜索空间，搜索效率更高；通过变形检验可以判断历元间的坐标差是由位移引起还是由误差引起，提高了处理结果的可靠性；预报精度更高，提高了GPS变形监测数据的处理效率。

附图说明

在下文中将基于仅为非限定性的实施例并参考附图来对本发明进行更详细的描述。其中：

图1为根据本发明的方法的流程框图。

具体实施方式

下面将参照附图来详细地介绍本发明。

本发明的目的在于克服现有技术中的用于处理GPS变形监测数据的方法的不足，提供一种基于变形检验的GPS单历元变形监测数据处理方法。

图1显示了根据本发明的方法的流程图。

参照图1，根据本发明的方法主要包括如下步骤：

第一步：获取基准站和监测站当前历元观测数据及广播星历数据，利用相位观测值组成双差观测方程。

第二步：根据变形监测对象的特性，设置两历元间变形的可发生的最大变形量。

第三步：采用Cholesky分解构建模糊度搜索空间。

其中：

$-\frac{{\mathrm{\δd}}_{\max }}{{\mathrm{\λl}}_{11}}\≤{\mathrm{\δ}}_{\▿\mathrm{\Δ}{N}_1}\≤\frac{{\mathrm{\δd}}_{\max }}{{\mathrm{\λl}}_{11}}$

$\frac{-\sqrt{{{\mathrm{\δd}}_{\max }}^2-B^2}-l_{21}\mathrm{\λ}{\mathrm{\δ}}_{\▿\mathrm{\Δ}{N}_1}}{{\mathrm{\λl}}_{22}}\≤{\mathrm{\δ}}_{\▿\mathrm{\Δ}{N}_2}\≤\frac{\sqrt{{{\mathrm{\δd}}_{\max }}^2-B^2}-l_{21}\mathrm{\λ}{\mathrm{\δ}}_{\▿\mathrm{\Δ}{N}_1}}{{\mathrm{\λl}}_{22}}$

$\frac{-\sqrt{{{\mathrm{\δd}}_{\max }}^2-B^2-C^2}-l_{31}\mathrm{\λ}{\mathrm{\δ}}_{\▿\mathrm{\Δ}{N}_1}-l_{32}\mathrm{\λ}{\mathrm{\δ}}_{\▿\mathrm{\Δ}{N}_2}}{\mathrm{\λ}{l}_{33}}\≤{\mathrm{\δ}}_{\▿\mathrm{\Δ}{N}_3}\≤\frac{\sqrt{{{\mathrm{\δd}}_{\max }}^2-B^2-C^2}-l_{31}\mathrm{\λ}{\mathrm{\δ}}_{\▿\mathrm{\Δ}{N}_1}-l_{32}\mathrm{\λ}{\mathrm{\δ}}_{\▿\mathrm{\Δ}{N}_2}}{{\mathrm{\λl}}_{33}}$

其中：

$B=l_{11}\·\mathrm{\λ}{\mathrm{\δ}}_{\▿\mathrm{\Δ}{N}_1,}$ $C=l_{21}\·\mathrm{\λ}{\mathrm{\δ}}_{\▿\mathrm{\Δ}{N}_1}+l_{22}\·\mathrm{\λ}{\mathrm{\δ}}_{\▿\mathrm{\Δ}{N}_2},$

δd_max为监测点的最大变形量，为双差整周模糊度中误差，λ为载波相位的波长，

$L^{-1}=\left[\begin{array}{ccc}{l}_{11}& & \\ l_{21}& l_{22}& \\ l_{31}& l_{32}& l_{33}\end{array}\right]$ 为求解条件方程组和误差方程组所组成的对称正定的方程组（法方程）矩阵AA^T的Cholesky分解的下三角矩阵。

第四步：利用最小二乘法搜索模糊度，其中方差比（Ratio值）需要大于3。

第五步：利用最小二乘法获取监测点的单历元整周模糊度固定解坐标X、Y、Z以及其协方差矩阵。

第六步：利用平均间隙法对当前历元与前一历元间的变形值进行变形检验。

检验的具体过程如下：

①由两历元观测值改正数（残差）获取的经验方差：

$S_0^2=\frac{{\left(V^T\mathrm{PV}\right)}_I+{\left(V^T\mathrm{PV}\right)}_{\mathrm{II}}}{f}$

其中V为观测值残差向量，P为观测值权阵，f为两历元自由度之和，即f＝n₁-3+n₂-3，其中n₁,n₂为两历元双差观测方程的个数。

②由两历元的坐标差（即所谓间隙）d_i(i＝1，2，3)构成单位权经验方差

${\stackrel{\‾}{S}}_0^2=\frac{d^T{P}_{d}d}{3}$

d＝X_II-X_I，表示监测点前后两历元的坐标之差向量，也可表示为：

$d=\left(\begin{array}{cc}-1& 1\end{array}\right)\left(\begin{array}{c}{X}_I\\ X_{\mathrm{II}}\end{array}\right)$

P_d为d的权阵：

$P_d=Q_d^{-1}$

此处Q_d由 $d=\left(\begin{array}{cc}-1& 1\end{array}\right)\left(\begin{array}{c}{X}_I\\ X_{\mathrm{II}}\end{array}\right)$ 根据协方差传播定律得到：

$Q_d=\left(\begin{array}{cc}-1& 1\end{array}\right)\left[\begin{array}{cc}{Q}_{X_I{X}_I}& 0\\ 0& Q_{X_{\mathrm{II}}{X}_{\mathrm{II}}}\end{array}\right]\left(\begin{array}{c}-1\\ 1\end{array}\right)=Q_{X_I{X}_I}+Q_{X_{\mathrm{II}}{X}_{\mathrm{II}}}$

③构建统计量：

$F_{h,f}=\frac{{\stackrel{\‾}{S}}_0^2}{S_0^2}=\frac{d^T{P}_{d}d}{h{S}_0^2}$

此统计量服从F分布，它的自由度分别为的自由度3和的自由度f。

④位移判断。选用显著水平α=0.01，进行右尾检验判断是否大于即将算出的F_h,f和从F分布表中查出的F_α(h,f)分位值比较，如果

F_h,f＞F_α(h,f)

则表明大于即有：

P{[F_h,f＞F_α(h,f)]}＝α

说明两历元间有位移发生，否则，表明两历元间没有发生位移。

第七步：经第六步的检验后两历元间未发生位移，则将两个历元的观测方程叠加的抗差最小二乘静态解作为当前历元的结果；

若第六步的检验结果为两历元间有位移发生，则通过抗差序贯平差获取当前历元结果。

据本发明的用于处理GPS变形监测数据的方法相对现有技术带来了如下进步：避免了周跳探测；采用Cholesky分解构建模糊度搜索空间，搜索效率更高；通过变形检验可以判断历元间的坐标差是由位移还是由误差引起，提高了处理结果的可靠性；预报精度更高。

虽然已经参考优选实施例对本发明进行了描述，但在不脱离本发明的范围的情况下，可以对其进行各种改进并且可以用等效物替换其中的部件。本发明并不局限于文中公开的特定实施例，而是包括落入权利要求的范围内的所有技术方案。

Claims (5)

1.一种用于处理GPS变形监测数据的方法，包括：

第一步：获取基准站和监测站当前历元观测数据及广播星历数据，利用相位观测值组成双差观测方程；

第二步：设置两历元间的最大变形量；

第三步：采用Cholesky分解构建模糊度搜索空间；

第四步：搜索模糊度，并且使方差比大于3；

第五步：获取监测点的单历元整周模糊度固定解坐标X、Y、Z和其协方差矩阵；

第六步：利用平均间隙法对当前历元和前一历元间的变形值进行变形检验；

第七步：经第六步的检验后，如果两历元间未发生位移，则将两个历元的双差观测方程叠加进行抗差最小二乘静态解算获得当前历元的结果；若第六步的检验结果为两历元间有位移发生，则根据抗差序贯平差获取当前历元的结果。

2.根据权利要求1所述的方法，其特征在于，在第三步中：

$-\frac{{\mathrm{\δd}}_{max}}{{\mathrm{\λl}}_{11}}\≤{\mathrm{\δ}}_{\▿{\mathrm{\ΔN}}_1}\≤\frac{{\mathrm{\δd}}_{max}}{{\mathrm{\λl}}_{11}}$

$\frac{-\sqrt{{{\mathrm{\δd}}_{\max }}^2-B^2}-l_{21}{\mathrm{\λ\δ}}_{\▿{\mathrm{\ΔN}}_1}}{{\mathrm{\λl}}_{22}}\≤{\mathrm{\δ}}_{\▿{\mathrm{\ΔN}}_2}\≤\frac{\sqrt{{{\mathrm{\δd}}_{\max }}^2-B^2}-l_{21}{\mathrm{\λ\δ}}_{\▿{\mathrm{\ΔN}}_1}}{{\mathrm{\λl}}_{22}}$

$\frac{-\sqrt{{{\mathrm{\δd}}_{\max }}^2-B^2-C^2}-l_{31}{\mathrm{\λ\δ}}_{\▿{\mathrm{\ΔN}}_1}-l_{32}{\mathrm{\λ\δ}}_{\▿{\mathrm{\ΔN}}_2}}{{\mathrm{\λl}}_{33}}\≤{\mathrm{\δ}}_{\▿{\mathrm{\ΔN}}_3}\≤\frac{\sqrt{{{\mathrm{\δd}}_{\max }}^2-B^2-C^2}-l_{31}{\mathrm{\λ\δ}}_{\▿{\mathrm{\ΔN}}_1}-l_{32}{\mathrm{\λ\δ}}_{\▿{\mathrm{\ΔN}}_2}}{{\mathrm{\λl}}_{33}}$

其中：

$B=l_{11}\·{\mathrm{\λ\δ}}_{\▿{\mathrm{\ΔN}}_1},C=l_{21}\·{\mathrm{\λ\δ}}_{\▿{\mathrm{\ΔN}}_1}+l_{22}\·{\mathrm{\λ\δ}}_{\▿{\mathrm{\ΔN}}_2},$

δd_max为监测点的最大变形量，δ_▽ΔN为双差整周模糊度中误差，λ为载波相位的波长，

$L^{-1}=\left[\begin{array}{ccc}{l}_{11}& & \\ 1_{21}& l_{22}& \\ l_{31}& 1_{32}& l_{33}\end{array}\right]$ 为求解条件方程组和误差方程组所组成的对称正定的方程组矩阵AA^T的Cholesky分解的下三角矩阵。

3.根据权利要求1或2所述的方法，其特征在于，在第六步中还包括如下子步骤：

第一子步骤：通过两历元观测值改正数获取经验方差：

$S_0^2=\frac{{\left(V^{T}PV\right)}_I+{\left(V^{T}PV\right)}_{II}}{f}$

其中f为两历元自由度之和，即f＝n₁-3+n₂-3，其中n₁,n₂为两个历元的双差观测方程的个数，V为观测值残差向量，P为观测值权阵；

第二子步骤：由两历元的坐标差d_i(i＝1，2，3)获取单位权经验方差 ${\stackrel{\‾}{S}}_0^2:$

${\stackrel{\‾}{S}}_0^2\frac{d^T{P}_{d}d}{3}$

其中，d＝X_Ⅱ-X_Ⅰ，表示监测点在前后两历元的坐标之差向量：

$d=\left(\begin{array}{cc}-1& 1\end{array}\right)\left(\begin{array}{c}{X}_I\\ X_{II}\end{array}\right)$

其中P_d为d的权阵：

$P_d=Q_d^{-1}$

此处Q_d由 $d=\left(\begin{array}{cc}-1& 1\end{array}\right)\left(\begin{array}{c}{X}_I\\ X_{II}\end{array}\right)$ 根据协方差传播定律得到：

$Q_d=\left(\begin{array}{cc}-1& 1\end{array}\right)\left[\begin{array}{cc}{Q}_{X_I{X}_I}& 0\\ 0& Q_{X_{II}{X}_{II}}\end{array}\right]\left(\begin{array}{c}-1\\ 1\end{array}\right)=Q_{X_I{X}_I}+Q_{X_{II}{X}_{II}}$

第三子步骤：构建统计量：

$F_{h,f}=\frac{{\stackrel{\‾}{S}}_0^2}{S_0^2}=\frac{d^T{P}_{d}d}{{\mathrm{hS}}_0^2}$

第四子步骤：选用显著水平α＝0.01，将获得的F_h,f与从F分布表中查出的F_α(h,f)分位值比较，如果

F_h,f＞F_α(h,f)

则表明大于即有：

P{[F_h,f＞F_α(h,f)]}＝α

从而判断出两历元间有位移发生；

如果无法得出F_h,f＞F_α(h,f)，则判断出两历元间没有发生位移。

4.根据权利要求1或2所述的方法，其特征在于，在第四步中，利用最小二乘法来搜索模糊度。

5.根据权利要求1或2所述的方法，其特征在于，在第五步中，利用最小二乘法来获取监测点的单历元整周模糊度固定解坐标X、Y、Z和其协方差矩阵。