DOI:10.19476/j.ysxb.1004.0609.2002.05.035

矿区开采沉陷GPS快速观测数据预处理方法

韩保民 欧吉坤 柴艳菊 卢秀山

  中国科学院测量与地球物理研究所动力大地测量学开放研究实验室  

  山东科技大学地球信息科学与工程学院 武汉430077  

  泰安  

摘 要：

随着精密GPS快速定位技术被广泛应用于矿区开采沉陷的日常观测 , 数据处理变得越来越重要 , 其质量好坏将直接影响到GPS快速定位精度 , 而粗差和周跳的探测与修复又是其中最重要组成部分。鉴于拟准检定法可以同时发现、定位测量数据中的多个粗差并估算其大小和精度 , 且予以修正 , 所以用粗差的拟准检定法来对GPS快速定位数据进行处理。首先概述了用拟准检定法处理GPS快速定位的数学模型 , 然后给出了简单的实施步骤。为了便于比较 , 还简单介绍了常用的抗差估计及统计检验的数学模型 , 最后通过 1个算例验证了本方法的可行性和有效性 , 并和另 2种方法作了比较

关键词：

GPS快速定位;粗差和周跳;数据处理;拟准检定法;

中图分类号： P228.4;TD327

收稿日期：2001-12-21

基金：国家自然科学基金资助项目 ( 4 0 0 740 0 3, 40 14 4 0 18);中科院知识创新工程项目 (KZCX2 -10 6 );

Method for processing data observed from GPS for subsidence surveying in mining area

Abstract：

The precision of GPS rapid positioning is directly related to the quality of data processing, so the data processing is an important procedure in the GPS rapid positioning. Since the "quasi-accurate detection of gross errors (QUAD) " method has the character of simultaneously detecting, locating multiple gross errors in surveying data, estimating their values and accuracies and then correcting them, the QUAD method is introduced to process GPS rapid positioning data. The mathematics model is summarized. Then detailed implementation procedure is presented. In order to make a comparison, the models of the other two methods like robust estimation and hypothesis test are also summarized. At last, one numerical example is given out to illustrate the feasibility and validity of this method.

Keyword：

GPS rapid positioning; gross error and cycle slip; data processing; QUAD method;

Received： 2001-12-21

近几十年来, 随着科学技术的迅猛发展, 特别是各种地面高精度观测技术及以全球定位系统 (GPS) 为代表的空间大地测量技术的突破性进展, 大大丰富了矿区开采沉陷的观测手段。 由于具有高效、快速、全天候、操作方便等特点, GPS快速定位已被广泛应用于开采沉陷的日常观测和对开采区内各种建筑物、构筑物的变形监测中。准确快速地求解模糊度是GPS快速、精密定位的关键性问题。 一旦整周模糊度确定, 相位观测值即转变为较精确的距离观测值, 由此可以求出精确的测点坐标 ^[1] 。但大多模糊度搜索方法易受粗差和周跳的影响, 因此, 数据预处理是GPS快速定位中重要一环, 其质量好坏将直接影响到GPS快速定位精度的高低, 而粗差和周跳的探测与修复又是其中最重要组成部分 ^[2,3,4,5,6] 。 对单频GPS快速定位来说, 如果仅用到L1频率的几个历元相位观测值, 这会给快速定位的数据质量控制带来较大的困难。如常用的双频P码或根据电离层延迟变化、多项式拟合法或高次差法和小波方法等来探测和修复GPS相位观测值中的粗差和周跳, 由于种种原因而变得很困难。 因而, 在实际GPS快速定位的数据处理中常用抗差估计及统计检验的方法 ^[7,8,9] 。 鉴于拟准检定法可以同时辨识、定位并估算观测数据中的多个粗差 ^[10,11,12] , 本文作者尝试将其用于探测和修复GPS相位观测值中的粗差和周跳。 首先分析了用这种方法进行数据处理的数学模型, 最后用一个算例和常用的抗差估计及统计检验等方法作了比较, 结果验证了本方法的可行性和有效性。

1 数学模型

1.1 拟准检定法的数学模型

设某历元两GPS接收机共视k+1颗GPS卫星, 则可组成k个双差观测方程, 其线性化后的双差观测方程为

$(\mathit{A}{}_i,\mathit{B})\left(\begin{array}{l}\stackrel{—}{{\displaystyle \mathit{X}}}\\ \mathit{{\rm N}}\end{array}\right)=\mathit{L}{}_i+\mathit{\Delta }{}_i\text{?}\text{?}\text{?}(1)$

式中 i为历元号; A_i为k×3维系数矩阵; $\stackrel{—}{{\displaystyle \mathit{X}}}$ 为未知基线向量改正值组成的向量, $\stackrel{—}{{\displaystyle \mathit{X}}}=(\delta X,\delta Y,\delta {\rm Z}){}^{\text{Τ}}$ ; N为整周模糊度向量; L_i为单频双差观测值向量; Δ_i为真误差向量; 观测值的权阵为P。

设这两台GPS接收机连续有n个历元共视了这k+1颗卫星, 由于各个历元的模糊度参数及其系数和未知基线向量改正参数不变, 则这n个历元总的双差观测方程为

可写成

AX=L+Δ (2b)

式中 A为n1×m维 (n1=k×n, m=k+3) 系数矩阵, X为m维待估参数向量。

记 R=I-A (A^TPA) ^-1A^TP, 容易导出如下确定关系式:RΔ=-RL。 由于R的秩为n-m, 是秩亏矩阵。 为此, 可选出r (>m) 个拟准观测, 然后附加“拟准观测的真误差范数极小”的条件, 即:

$\parallel \hat{\Delta }{}_r\parallel {}_{{\rm P}}^2=\Delta {}_r^{\text{Τ}}\mathit{{\rm P}}\Delta {}_r=\min \text{?}\text{?}\text{?}(3)$

此时取: $\mathit{G}=\mathit{G}{}_Q=(0,\mathit{A}{}_r^{\text{Τ}})?\mathit{{\rm P}}=\mathit{{\rm P}}{}_Q+\left(\begin{array}{cc}0& 0\\ 0& \mathit{{\rm P}}{}_r\end{array}\right)$ , 则附加条件可表示为

G_QP_QΔ_Q=Α^T_rP_rΔ_r=0 (4)

可由方程RΔ=-RL求解真误差估值 $\hat{\Delta }$ :

$\hat{\Delta }{}_Q=\left(\begin{array}{l}\hat{\Delta }{}_l\\ \hat{\Delta }{}_r\end{array}\right)=-(\mathit{R}{}^{\text{Τ}}\mathit{R}+\mathit{{\rm P}}{}_Q\mathit{G}{}_Q^{\text{Τ}}\mathit{G}{}_Q\mathit{{\rm P}}{}_Q){}^{-1}\mathit{R}{}^{\text{Τ}}\mathit{R}\mathit{L}\text{?}\text{?}\text{?}(5)$

通过初选、复选拟准观测求得 $\hat{\Delta }?$ 根据 $\hat{\Delta }$ 呈现明显分群的特征, 可以确定粗差或周跳的位置。 最后用如下公式来估计粗差或周跳的大小:

$\widehat{?}{}_Q$ = (C^T_bPRC_b) ^-1C^T_bPRL (6)

1.2 双因子抗差估计模型

采用杨元喜等提出的双因子抗差估计模型 ^[8,13] , 构造的相关等价权元素为

$\begin{array}{l}\stackrel{—}{{\displaystyle {\rm P}}}{}_{ij}={\rm P}{}_{ij}\gamma {}_{ij},\gamma {}_{ij}=\sqrt{\gamma {}_{ii}\gamma {}_{jj}}?\text{?}\text{?}\text{?}(7)\\ \gamma {}_{ii}=\{\begin{array}{l}\text{?}\text{?}\text{?}\text{?}1\text{?}\text{?}\text{?}\text{?},|\stackrel{?}{{\displaystyle \upsilon }}{}_i|\le k{}_0\\ \frac{k{}_0}{|\stackrel{?}{{\displaystyle \upsilon }}{}_i|}\left(\frac{k{}_1-|\stackrel{?}{{\displaystyle \upsilon }}{}_i|}{k{}_1-k{}_0}\right){}^2\text{?},k{}_0＜|\stackrel{?}{{\displaystyle \upsilon }}{}_i|\le k{}_1\\ \text{?}\text{?}\text{?}\text{?}0\text{?}\text{?}\text{?}\text{?},|\stackrel{?}{{\displaystyle \upsilon }}{}_i|＞k{}_1\end{array}\text{?}\text{?}\text{?}(8)\end{array}$

式中 P_ij为观测值的权, γ_ii和γ_jj为等价权因子, $\stackrel{?}{{\displaystyle \upsilon }}{}_i$ 为标准化残差, k₀和k₁是常数。 当观测误差落入正常段时, $\stackrel{—}{{\displaystyle {\rm P}}}{}_{ij}={\rm P}{}_{ij}$ (等权) ; 当观测误差落入可疑段又未呈现显著异常时, $\stackrel{—}{{\displaystyle {\rm P}}}{}_{ij}＜{\rm P}{}_{ij}$ (降权) ; 当观测误差呈现显著异常时, $\stackrel{—}{{\displaystyle {\rm P}}}{}_{ij}=0$ (剔除) 。

1.3 假设检验模型

采用B检验 ^[3] , 其统计量为

T_i $=(\mathit{{\rm P}}\mathit{V}){}_i/(\sigma {}_0\sqrt{(\mathit{{\rm P}}\mathit{R}){}_{ii}})\text{?}\text{?}\text{?}(9)$

式中 R为平差因子阵, P为权阵, σ₀为单位权中误差。 如果最大统计量T_i>K_α=F ${}_{\alpha }^{1/2}$ =F ${}_{\alpha }^{1/2}$ (1, ∞, 0) , 则认为该观测值有粗差或周跳, 在剔除或修正该观测后重新进行统计检验, 直到所有统计量都通过检验为止。 本文中取置信水平α=0.001, 阀值K_α=3.29 (常用) 。

2 算例及分析

试验描述:在某矿一条长为4.7 km的岩移观测线上进行实验观测, 采样间隔为15 s, 截止高度角为15°, 共观测了5个历元。 组成双差时的5个卫星对为:11-29, 18-29, 19-29, 22-29和28-29。 为了更好的说明问题, 我们在第一历元的28-29卫星对、第三历元的22-29卫星及第五历元的28-29卫星对上分别模拟了一个0.12 m, 0.15 m和0.22 m的粗差, 在第4历元的11-29, 18-29卫星对上分别模拟了1周的周跳。 由于周跳具有继承性, 因此在下一历元的相应卫星对上也有1周的周跳 (模拟值见表1的4列) 。 下面结合算例来简要说明用拟准检定法探测和修复GPS相位观测值中粗差和周跳的步骤。

先求出残差V=-RL, 然后经过一次初选和数次复选拟准观测, 根据真误差估值存在的明显分群现象来判定那些含粗差或周跳的观测值, 并给出较为准确的估值。 在此过程中, 先要计算初选拟准观测时用到的分类指标u_i, 然后还要计算各观测值的真误差估值 $\hat{\Delta }{}_i$ 及复选拟准观测用到的分类指标W^Ⅰ_i和W^Ⅱ_i。 其中, $W{}_i^{\mathrm{Ⅰ}}=\hat{\Delta }{}_i/C{}_1?C{}_1=1.483\underset{i}{{\displaystyle \text{m}\text{e}\text{d}}}|\hat{\Delta }{}_i|?(i=1,?,n)?W{}_i^{\mathrm{Ⅱ}}=\hat{\Delta }{}_i/C{}_3?C{}_3=\sqrt{\hat{\Delta }{}_r^{\text{Τ}}\hat{\Delta }{}_r/(r-1)},r$ 为拟准观测数。

从最后求得的 $\hat{\Delta }{}_i?|W{}_i^{\mathrm{Ⅰ}}|$ 和|W^Ⅱ_i|可知, 后面7个真误差估值及对应的|W^Ⅰ_i|和|W^Ⅱ_i|远远大于其它值, 存在明显分群现象, 可以认为它们所对应的观测值确实存在粗差或周跳。 从按式 (6) 计算的粗差估值来看, 也和模拟值很接近。 最后只要从观测值中减去这些估值就达到修复周跳的目的。 至此, 本算例模拟的粗差和周跳都被探测出来, 并得到较好的修正。 有关计算结果见表1、图1。 表1中列1为历元号, L₀为没有粗差和周跳的观测值, ?为模拟的粗差和周跳, L=L₀+?, 表示含粗差或周跳的双差观测值。 从图中可更直接地看出:经过初选和2次复选拟准观测, 真误差估值分群的现象表现得更突出, 含粗差或周跳的真误差估值明显离群漂浮到上面, 而其他没有“问题”的真误差估值则沉到底下。

为了便于比较, 文中同时用抗差估计模型和统计检验方法对上述数据进行处理, 其各项指标比较结果见表2, 处理后的残差图见图2。 从表2、图2中不难看出, 拟准检定法在对污染率较高的GPS快速定位数据 (28%) 的处理时明显优于其他2种方法, 不仅定位准确, 而且估值也很接近模拟值, 经改正后的观测值的残差也趋近0。

表1 用拟准检定法对快速定位的数据预处理结果

Table 1 Processing results obtained from using QUAD method for pre-processing rapid positioning data

1	2	3	4	5	6	7	8	9	10	11	12	13	14	15
Epoch	No.	L0/m	?/m	L/m	No.	Class	ui	No.	$$|\hat{\Delta }{}^{(1)}|$$	|WⅠ (1) |	No.	$$|\hat{\Delta }{}^{(2)}|$$	|WⅠ (2) |	|WⅡ (1) |
1	1	-1.155	0.000	-1.155	13	2	0.015	19	0.000	0.000	15	0.000	0.065	0.103
	2	-1.780	0.000	-1.780	18	2	-0.280	9	0.000	0.001	12	0.000	0.069	0.109
	3	-2.319	0.000	-2.319	8	2	0.308	6	0.000	0.001	24	0.001	0.076	0.119
	4	-2.911	0.000	-2.911	6	2	0.445	10	0.000	0.001	4	0.001	0.112	0.177
	5	-3.478	0.120	-3.358	9	2	0.474	7	0.000	0.001	13	0.001	0.180	0.284
	6	-1.164	0.000	-1.164	19	2	0.527	18	0.000	0.001	23	0.001	0.181	0.286
	7	-1.788	0.000	-1.788	7	2	0.593	8	0.000	0.002	11	0.001	0.218	0.344
	8	-2.328	0.000	-2.328	10	2	0.769	13	0.000	0.002	3	0.002	0.281	0.443
	9	-2.917	0.000	-2.917	15	2	1.014	24	0.006	0.166	2	0.002	0.312	0.493
2	10	-3.485	0.000	-3.485	16	2	-1.040	23	0.010	0.259	19	0.002	0.339	0.535
	11	-1.159	0.000	-1.159	17	2	-1.170	4	0.011	0.291	1	0.003	0.401	0.633
	12	-1.781	0.000	-1.781	12	2	1.175	3	0.015	0.393	18	0.004	0.569	0.898
	13	-2.322	0.000	-2.322	11	2	1.184	15	0.026	0.674	20	0.005	0.674	1.065
	14	-2.911	0.150	-2.761	20	2	1.364	12	0.026	0.680	9	0.006	0.825	1.302
3	15	-3.474	0.000	-3.474	25	2	-1.596	16	0.036	0.955	6	0.006	0.920	1.452
	16	-1.154	0.190	-0.964	5	2	-2.689	2	0.041	1.080	7	0.007	0.980	1.548
	17	-1.775	0.190	-1.585	14	2	-2.209	21	0.050	1.325	8	0.007	1.005	1.587
	18	-2.317	0.000	-2.317	1	1	-0.423	20	0.057	1.491	10	0.010	1.397	2.206
	19	-2.907	0.000	-2.907	2	1	-0.994	11	0.077	2.030	5	0.118	17.148	27.078
4	20	-3.470	0.000	-3.470	3	1	1.480	1	0.091	2.391	14	0.149	21.789	34.407
	21	-1.159	0.190	-0.969	4	1	1.279	22	0.102	2.665	21	0.189	27.513	43.446
	22	-1.779	0.190	-1.589	21	1	-0.354	25	0.123	3.222	22	0.190	27.733	43.792
	23	-2.322	0.000	-2.322	22	1	-0.883	17	0.138	3.626	16	0.194	28.277	44.652
	24	-2.909	0.000	-2.909	23	1	-1.481	14	0.151	3.953	17	0.195	28.467	44.952
5	25	-3.475	0.220	-3.255	24	1	-1.274	5	0.163	4.281	25	0.219	32.012	50.549
							r=8		r=12			r=18

表2 用3种不同方法对快速定位数据处理后的各项指标比较

Table 2 Comparison of several index obtained from using three different methods for processing rapid positioning data

Method	Gross errors and cycle slips in raw data	Detected gross errors and cycle slips	Estimators of gross errors /m	Success rate	Rms./m	Iteration number
QUAD	5, 14, 16, 17, 21, 22, 25	5, 14, 16, 17, 21, 22, 25	0.118, 0.149 0.189 0.190 0.194 0.195 0.219	100%	0.0042	3
Robust estimation	See above	5, 14, 17, 25	Lower weights for gross errors	57.1%	0.548	26
Hypothesis test	See above	5, 14, 25	Lower weights for gross errors	42.9%	0.602	3

图1 拟准检定法处理时真误差值的分群现象

Fig.1 ‘Hive off’ phenomenon of real error estimators at processing using QUAD method

图2 3种不同方法处理后的残差示意图

Fig.2 Illustration of residuals obtained from three different methods

(△—Robust estimation residual; ●—QUAD residual; ○— Hypothesis test residual)

3 结论

结果表明, 拟准检定法思路简单、直观, 在对GPS快速定位数据作预处理时有以下特点:①适合于GPS单频接收机; ②可以同时探测和修复GPS相位观测值中的粗差和周跳, 并能给出比较准确的估值然后对其进行修正; ③对观测数据污染较重的情况也能作较好的处理。

参考文献

[1]　HofmannWellenhofB , LichteneggerH , CollinsJ .ClobalPositioningSystem, TheoryandPractice (Fourthrevisededition) [M ].NewYork:SpringerVerlag, 1994, 181-200.

[2]　SalzmannM .Areal timequalitycontrolprocedureforuseinintegratednavigationsystems[J].TheHydrographicJournal, 1995, 72:25-30.

[3]　FreiE , BeutlerG .Rapidstaticpostitionbassedonthefastambiguityresolutionapproach“FAFA”theoryandfirstresults[J], ManuscriptaGeodactica, 1990, 15:83-55.

[4]　MertikesSP , RizosC .On linedetectionofabruptchangesinthecarrier phasemeasurementsofGPS[J].JournalofGeodesy, 1997, 71:469-482.

[5]　TeurnssenPJG , KleusbergE .QualityControlandGPS , InGPSforGeodesy (2ndEdition) [M ].Berlin:SpringerVerlag, 1998, 271-316.

[6]　BisnathSB , KimD , LangleyRB .Anewapproachtoanoldproblem:carrier phasecycleslips[J].GPSWorld, 2001, 12 (5) :46-51.

[7]　LIDe ren.ErrorProcessingandReliabilityTheory[M ].Beijing:SurveyingandMappingPress, 1988.187-201.

[8] 　YANGYuan xi, SONGLi jie, XUTian he.Robustesti mationforcorrelatedGPSbaselinenetwork[A].Trans actionsonGeodesy[C].Beijing:SurveyingandMappingPress, 1999.32-39.

[9]　RemondiBW .UsingtheGlobalPositioningSystem (GPS) PhaseObservableforRelativeGeodesy:Modeling, ProcessingandResults[D].Austin:UniversityofTexas, 1984.

[10]　OUJi kun.Anewmethodofindentifyingandlocatinggrosserrors quasi accuratedetection[J].ScienceBul letin, 1999, 44 (16) :1777-1781.

[11]　OUJi kun.Quasi accuratedetectionofgrosserrors (QUAD) [J].ActaGeodaeticaetCartographicaSinica, 1999, 28 (1) :15-20.

[12] 　OuJi kun.Selectionofquasi accurateobservationand“hiveoff”phenomenaabouttheestimatorsofrealerrors[J].ActaGeodaeticaetCartographicaSinica.2000, 29 (Sup) :5-10.

[13]　CasparyW , BoruttaH .Robustestimationindeforma tionmodels[J].SurveyReview, 1987, 28 (223) :29-45.