基于QOCA软件的GNSS连续站中非构造噪声剔除方法研究

纪海源,何远梅,郝 明

(1.陕西工业职业技术学院,陕西 咸阳 712000; 2.中国有色金属工业西安勘察设计研究院有限公司,西安 710000; 3.中国地震局第二监测中心,西安 710043)

0 引言

连续GNSS变形监测是地壳形变监测的主要方法。研究发现,在连续站位移时间序列中存在着非地壳内部动力因素造成的非构造噪声,这种非构造噪声由很多因素引起,包括数据处理软件中模型没有消除的、地表质量负荷、随着季节性下雨与下雪等负荷及观测墩的季节性热胀冷缩等。姚宜斌等[1]指出,GLOBK 软件假定测站运动呈线性变化,与现实中观测到的测站运动不符,可能会引入粗差。QOCA软件在数据处理中能够兼顾测站运动的非线性变化,运用QOCA软件剔除GNSS连续站位移序列中的非构造噪声,可为地壳的形变研究提供重要的基础数据。QOCA(Quasi-Observation Combination)是空间大地测量资料处理软件,可对各种松弛约束的大地测站坐标及速度解(如准观测值)进行联合解算,获取地壳形变信息,不仅能够对全球定位系统及卫星激光测距等空间大地拟观测值进行解算,还可以与陆地大地测量拟观测值(电子激光测距、三角测量、水准测量等)进行联合解算,能够联合合成孔径雷达数据(SAR)、重力、地震活动及地表螺变等数据进行解算。通过计算获得测站坐标、测站速度、网参数及同震与震后形变参数,探测出包含在拟观测值中的异常值,采取稳健性分析直接处理噪声时间序列。

1 非构造噪声

GAMIT/GLOBK处理得到的GNSS连续站位移时间序列是相对精确的,GAMIT中所加载的各种改正模型在很大程度上降低了噪声的影响,但是不能完全消除,故这种位移序列结果中除了包含真实的地壳构造形变外,还包含了一些干扰噪声,图1为连续站GSJN的位移时间序列,在NEU三个方向明显能看见周期性成分,这不是由地壳内部动力因素造成的。在GNSS监测过程中产生的干扰统称为非构造噪声,这种非构造噪声由很多因素引起,包括GAMIT中模型没有消除的地表质量负荷,随着季节性下雨、下雪的负荷,观测墩的季节性热胀冷缩等[1]。地壳的运动不可能呈周期性变化(一般为线性),故去除其中的非构造噪声信息对于有效运用GNSS数据研究构造形变场来说至关重要。

图1 位移时间序列Fig.1 Displacement time series

2 analyze_tseri与pca模块的数据分析原理

analyze_tseri与pca是QOCA软件的主要组成模块,对非构造噪声的分析有重要的作用。

analyze_tseri为时间序列分析(time series analysis),是一种动态数据处理及统计方法,广泛应用于多领域,可将来自地表与地下的策动源产生的确定性地形变模式(pattern)以参数化方式估计出来。把时间序列中未知的信号看作一随机过程,利用自回归分析、最大熵谱估计、ARMA模型拟合等统计方法把隐藏的信号分离出来,通常将地形变时间序列用数学模型表示为:

x(t)=x0+v·(t-t0)+∑[Sisin(ωit)

pca(principal component analysis)称为主分量分析、经验正交函数分析(EOF)、特征向量分析,能够把随时间变化的台站网时间序列分解成时间域的主分量及空间域的特征向量,这种分解按照每个主分量贡献的功率(能量)来排列,一般前几个主分量的变化特征可以将区域性时间变化特征最大限度地表现出来,这几个主分量对应的特征向量反映了这些时间变化强弱的空间分布[2]。如果把台站网中每个台站的时间序列排列起来,每个台站时间序列为一列,可表示为如下m×n(通常情况下,m>n)的数字矩阵:

可以由奇异值分析分解成如下形式:X=UΠV,式中,U为m×m正交归一矩阵,为m×n准对角矩阵,V为n×n正交归一矩阵。地形变分析中最常见的情况是m>n,且X的秩为n。这时X的方差矩阵可表为:

3 QOCA软件剔除非构造噪声

为了剔除GNSS位移时间序列中的非构造噪声,必须获取这种噪声的时间序列,分析其中包含的信息。根据黄立人等[3]的研究资料,把噪声按照功率谱密度P(f)与其对应的噪声频率f之间的关系P(f)∝(f)k分类为:若k=0则对应的是白噪声,与频率没有多大关系,各种频率下噪声大小一样。若k=1,则对应的是闪烁噪声,功率谱的密度与频率成反比,有一定的波动。若k=2,则对应的噪声称为随机漫步噪声,有一定的随机性,波动性大。噪声性质的确定关键是获得噪声序列并求出它的功率谱密度(简称功率谱)与谱指数。对于连续的位移时间序列分量来说,可采用模型消减这种噪声影响[3]。为了得到观测噪声的时间序列,对位移分量每日解观测序列建立参数模型:

式中,ti(i=1…N)为以年为单位的时间,在这里以每日的单日解构成时间序列。a为初始位置,b为速率,s1,s2,s3、s4分别为年、半年周期项系数,g为偏移,h为震后速率变化,k为震后速率衰减(指数模型),H为阶梯(heavisidestep)函数,cme为共模误差,υ为误差。以上数学模型可通过QOCA软件的模块来实现。

3.1 analyze_tseri分离线性成分

使用analyze_tseri模块对GLOBK处理结果的36个连续站进行位移时间序列分析,拟合其中的线性项、周年、半周年项等参数,得到拟合残差。

analyze_tseri输入文件准备。jz_unf.drv:驱动文件,这是analyze_tseri模块最重要的一个文件,包含所有输入输出文件的路径、名字及各种选择项的参数。analyze_tseri会根据驱动文件中的参数来运行。文件由用户建立,采用自由格式,主要文件包括coord.soln(台站先验坐标和速度文件)、stxyz.list(输入数据文件的列表文件)、jzsite.list(台站文件)、jz_tseri.para(估计参数文件)。其中,jz_tseri.para文件给出了每个台站要估计的参数及每个参数的情况,这是程序最关键的文件,参数拟合与估计的好坏全看它给的是否合理。

运行analyze_tseri指令,给出驱动文件名字(analyze_tseri analyze_tseri.drv)。

结果文件。jz_unf.out:输出结果文件。头文件中必须给出输出文件名及路径。这个文件包含了所有估计参数的结果,其中以大写字母开头的一般是估计参数解。小写字母行的意义如下:lsq_velo(速度项解)、lsq_coor(台站坐标解)、coor_xyz与velo_xyz(坐标与速度估计在xyz坐标系中的解)、jz_unf.resi(残差文件)。jz_unf.resi给出了第一个残差文件的名字及路径,考虑为后面的区域滤波做准备,必须输出残差文件。图2是根据参差文件jz_unf.resi画的残差时间序列。

3.2 pca分离共模误差

由analyze_tseri处理得到的残差图2可以看出,analyze_tseri所得的残差时间序列中仍可见周期性成分及不规则波动。根据Wdowinski S等[4]的研究表明,残差中包含的这类非构造噪声称为共模误差(common modeerror,即cme),这类噪声存在区域相关性,来源尚不清楚,有可能来自卫星轨道误差、水体及大气质量负荷、参考框架定义的不确定性等[5],此外还存在单站误差,源于天线墩不稳定、多路径效应等[6]。Dong等[5]采用pca方法对台站网中是否存在cme制定了标准,如果大于50%的站出现明显的空间响应,且这种模式特征值与特征值总和大于1%,按照这个标准一般构成cme的主体部分是排在前面的几个主成分。SCIGN通过pca分析的残差坐标序列文件表明,cme具有空间特征,以同一方式在一定的区域内影响各站点。因为这种影响不是地壳构造活动造成的,故有必要从坐标时间序列中剔除[7-9]。采用pca模块对analyze_tseri的残差文件得到的残差坐标序列进行主成分分析。

pca分析文件准备。pca_jz.drv:pca的驱动文件,装载所有输入输出文件的路径及名字及各种选择项的参数,这是一个必须由用户建立的文件,有3种类型,即regional filtering、tectonic signal、seasonal signal,pca的内插方式是不一样的。对于regional filtering工作类型,程序要求内插尽量不影响由没有空缺观测的台站给出的前1～2个主分量,对于后面两种工作类型,程序要求内插尽量不影响由没有空缺观测的台站给出的主要构造运动模式即季节性变化模式。处理选择regional filtering,需准备台站先验坐标及速度文件coord.soln,台站文件site_list,输入数据文件的列表文件pca_jz.list。

运行pca指令,给出驱动文件名字pca pca_jz.drv即可。

输出结果如下:

pca_jz.out:输出结果文件,给出输出文件名及路径,主要储存各个分量分解的本征值及其功率。

pca_jz.cpt:主分量文件,给出主分量文件的名字及路径。输出12个主分量(单位a,m)按东西、南北、垂直向依次排列如下:

* time cpt1 cpt2 cpt3 cpt4 cpt5 cpt6 cpt7 cpt8 cpt9 cpt10 cpt11 cpt12

2011.00100000 -0.004025 -0.001041 -0.000234 0.001173 0.001535 0.000044 -0.005732 -0.002654 -0.000331 -0.000613

-0.002197 0.000410

pca_jz.seign:输出特征向量文件,给出空间响应本征矢文件名字及路径。输出的前10个特征向量(无单位)按东西、南北、垂直向依次排列:

* long lati site cpt1 cpt2 cpt3 cpt4 cpt5 cpt6 cpt7 cpt8 cpt9 cpt10

104.6046 35.5543 GSDX_GNSS 0.758543 0.201866 0.131942 -0.020168 0.021891 -0.002359 -0.057869 -0.030266 0.434962 0.214931

pca_jz.res:滤波后残差序列输出文件,必须给出输出文件名及路径。通过pca的主分量分析得到了两个重要文件,即pca_jz.cpt、pca_jz.seign,包含主要的共模误差,可作为analyze_tseri的输入文件。

3.3 analyze_tseri剔除非构造噪声

经过pca主成分分析得到了主分量文件pca_jz.cpt及本征矢文件pca_jz.seign,将其加入到analyze_tseri的驱动文件中,可以剔除pca分离出来的共模误差。对剔除共模误差后的位移时间序列中的周年项、半周年项等非构造噪声参数信息进行拟合,从位移序列中剔除。QOCA会根据剔除干净的位移时间序列来拟合一个最终速度,这个速度即GNSS连续站在扣除非构造噪声后相对真实的线性运动速度。analyze_tseri处理应注意在驱动文件jz_flt.drv中的关键字cme_correction file行后面添加滤波文件pca_jz.seign pca_jz.cpt。

3.4 QOCA处理结果分析

经过QOCA软件对GNSS连续站的处理,已经将大量的非构造噪声剔除,得到需要的各种估计参数及速度文件。表1是从最终的结果文件jz_flt.out中提取的部分连续站的速度。

表1 剔除非构造噪声后的部分连续站速度统计Tab.1 Partial continuous station velocity statistics after removing unstructured nois

图3是用连续站GSJN剔除线性成分与共模误差后的残差文件画的,从图中明显能看到比滤波前的波动小,但是出现个别点误差的增大,这可能是因为个别单日解的误差较大。通过pca方法得到的结果与Dong等[7]的结果基本一致。

图3 剔除线性成分和共模误差后的GSJN站残差Fig.3 Residual of GSJN station after removing linear component and common mode error

将剔除非构造噪声前后的速度进行对比,表2列出了部分连续站的数据,在计算速度差均方根时采用全部36个测站,经过对比发现,在南北方向上最大速度差为0.6 mm,最小为0 mm,速度差均方根为0.26 mm。在东西方向上最大差0.5 mm,最小差0 mm,速度差均方根0.19 mm。在垂直方向上最大-1.8 mm,最小为-0.3 mm,速度差均方根0.64 mm。水平方向上的速度差均方根为0.33 mm,可见,垂直方向上速度差均方根是水平方向上的1.97倍。这符合GNSS高程测量精度小于水平方向上的规律。QOCA处理结果在水平方向还不是太明显,而在垂直方向上对提高GNSS高程测量有很大的作用。

表2 剔除非构造噪声前后部分连续站速度对比Tab.2 Velocity comparison of partial continuous stations before and after removing unstructured noise

4 结论

使用QOCA软件的analyze_tseri与pca模块对36个GNSS连续站位移时间序列进行分析处理,剔除连续站中的共模误差,获取了位移时间序列中的周期性参数信息,将其中的周期性等非构造噪声剔除,令GNSS连续站水平运动速度的精度提高了约0.3 mm/yr,垂直方向运动速度的精度提高了约0.6 mm/yr,得到相对真实的地壳形变信息,为地壳的形变研究提供了重要的基础数据。