GAMIT与GLOBK软件介绍

1.GAMIT软件功能及组成

(1)GAMIT软件模块

GAMIT软件由许多功能不同的模块组成，这些模块可以独立运行。各个模块具有一定的独立性，但它们之间又紧密地联系在一起，共同完成数据处理和分析的全过程。这些模块可按其功能分成两个部分：数据准备和数据处理。此外，该软件还带有功能强大的SHELL程序。数据准备部分包括原始观测数据的格式转换、计算卫星和接收机钟差、星历的格式转换等；数据处理部分包括观测方程的形成、轨道的积分、周跳的修复和参数的解算等。

1)数据准备模块

①MAKEXP：数据准备部分的驱动程序，建立所有准备文件的输出及一些模块的输入文件；

②BCTOT(NGSTOT)：将星历格式(RINEX、SP3、SP1)转换成GAMIT所需的卫星星历T文件；

③MAKEJ：读取观测数据，生成卫星钟差文件J文件；

④MAKEX：将原始观测数据的格式(RINEX)转换成GAMIT所需的接收机时钟文件K文件和观测文件X文件。

2)数据处理模块

①FIXDRV：数据处理部分的驱动程序。

②ARC：轨道积分模块；

③MODEL：求偏导数，组成观测方程；

④AUTCLN：进行相位观测值周跳和粗差的自动探测与修复；

⑤SINCLN：单站自动修复周跳；

⑥DBLCLN：双差自动修复周跳；

⑦CVIEW：在可视化界面下人工交互式修复周跳；

⑧CFMRG：用于创建SOLVE模块所需的M文件，选择和定义有关参数；

⑨SOLVE：利用双差观测值，按最小二乘法解算参数。

3)辅助模块

辅助模块包括CTOX、XTORX、TFORM等。

(2)GAMIT常用文件及格式说明

1)测站信息文件

接收机和天线的型号、版本、天线高等情况均记录于测站信息文件station.info中，该文件会被MAKEXP、MAKEX、MODEL模块读取。此文件由用户自己准备，当然也可自动生成，具体形式见表3-5。

表3-5　测站信息文件

表中上下两行内的各项相互对应，其中SITE为四个字符的测站名，Station Name为该测站的完整测站名，Session Start与Session Stop标明了该测站观测时段的起始与终止时间，Ant Ht为测站天线高，HtCod为测站天线高的量测类型，Ant N为天线北方向改正，Ant E为天线东方向改正，Receiver Type为接收机型号，Antenna Type为接收机天线型号，Dome是指是否有天线罩。测站天线高和天线高量测方式的输入需认真核对，一旦出错，将会使解算结果在高程方向上产生系统性偏差。

2)测站信息控制文件

测站信息控制文件sittbl.分为long型和short型两种文件格式。short型sittbl.文件适用于对GAMIT软件尚不够了解的初学者，该文件将大部分参数设定为缺省设置；而long型的sittbl.文件其可操作性更强，更有利于那些对GAMIT数据处理有一定了解的使用者根据具体情况修改相应参数。本文现仅以long型的sittbl.文件进行说明，文件格式见表3-6。

表3-6　测站信息控制文件

测站信息控制文件sittbl.中，SITE为四字符点名；FIX决定该测站是否为固定点，YYY表示是固定点，NNN表示不是固定点；WFILE表示是否存在水汽辐射计文件，若存在则为WVR文件名，一般情况设为NONE；COORD.CONSTR.表示测站三维坐标约束量，如表3-6中的0.005 0.005 0.010分别表示WUHN站的三维坐标约束量，单位：m。因为该例中将WUHN站作为固定点，即起算点，所以约束量较小，如果是非固定点的测站，约束量通常取20.00 20.00 20.00；EPOCH指参加计算的起始历元数，001-∗表示所有历元；CUTOFF指截止高度角，通常取15.0，单位：°；CLK指是否解算接收机钟差的漂移量；KLOCK、CLKFT指接收机钟差改正模型；DZEN为对流层干项延迟的计算模型，默认为SAAS；WZEN为对流层湿项延迟的计算模型，默认为SAAS；DMAP与NMAP分别指干项延迟映射因子和湿项延迟映射因子；MET.VALUE为标准气象参数，即气压、温度以及相对湿度。对于未在sittbl.文件中进行设置的测站，其测站参数将被赋予sittbl.文件中所指定的默认值。

3)测段信息控制文件

测段信息控制文件sestbl.主要是对GAMIT软件进行参数设定，其格式见表3-7(只列出主要部分)。

表3-7　测段信息控制文件

该文件中，Type of analysis指对解算方法进行选择，具体如下：

0-ITER：ARC(optional)，MODEL，AUTCLN(optional)，SOLVE，

SCANDD；

1-ITER：指两个0-ITER序列，但第一个是QUICK解；

2-ITER：指在1-ITER基础上再加一个序列，用于确定轨道。

Choice of observable指对观测量类型进行选择，具体如下：

LC_HELP：用LC观测值解算，利用虚拟电离层观测值解模糊度；

L1_ONLY：仅仅使用L1观测值，对于几公里的小网；

L2_ONLY：仅仅使用L2观测值，对于几公里的小网；

L1，L2_INDEPENDENT：联合使用L1、L2进行解算，对于几公里的小网。

Choice of experiment指对解算类型进行选择，具体如下：

RELAX：包括定位、定轨，并解地球自转参数(ERP)；

BASELINE：仅仅是定位。

Zenith Delay Estimation是指是否估算对流层天顶延迟，若选“YES”，则对下列项进行设置：

Interval zen：设置间隔多少小时估计一个天顶延迟参数；

Zenith Model：天顶延迟估算模型选择，PWL指“线形插值法”；

Atmospheric gradients：设置是否估算大气的水平梯度，默认为“NO”。一般来说，应该估算对流层天顶延迟，并按2小时估计一个天顶延迟参数。

4)坐标文件

GAMIT的输入L文件是站坐标文件，包括测站的先验坐标(近似坐标)，目前测站坐标以空间直角坐标表示。利用程序glbtol将apr文件转换为当前观测历元的L文件，而IGS站的apr文件可由ITRF直接获得。如果不利用GLOBK进行平差，可以直接将更新后的L文件作为当天的坐标平差结果。

测站近似坐标的正确与否对于基线解算精度有着较大的影响。在批处理中，其近似坐标是根据所读取的测站观测O文件自动生成。将高精度已知点或与其进行了长时间联测的点位坐标作为基线处理的参考基准进行约束，解算各站点间的基线结果。在此必须强调，近似坐标文件所提供的各站点近似坐标其绝对误差必须小于300m。这是因为在每个历元相位观测值的处理中我们都需要计算出接收机的钟差。如果接收机和卫星的坐标已知，根据所观测的伪距值，可按如下公式求出接收机钟的偏差：

式中，p1为伪距观测值，r为卫星与接收机间的预报距离，c为光速。当接收机钟的精度达到1微秒时，据此所解得的基线解算精度约为1mm。为了达到这样一个精度水平，我们要求测站近似坐标误差对接收机钟所造成的偏差影响要小于1微秒。此时，则要求测站的近似坐标误差不能大于300m。

5)基线处理结果文件和站坐标系

GAMIT的数据处理输出文件为基线约束解O文件和基线松弛解h文件，主要包括测站的球面坐标及基线结果。球坐标与ITRF坐标的转换公式为：

使用球坐标系的优点是：球坐标系与直角坐标系之间的转换简单，不需要像大地坐标那样进行迭代运算，经度与大地经度相同，纬度与大地纬度比较接近，并且径向方向的变化可以近似认为是高程的变化。严格地说，GNSS数据处理是为了解算未知点的三维坐标，而不是基线分量。基线是由坐标计算的，即先有坐标后有基线，一些测量工作者往往混淆了这一概念。在建立误差方程时，一般是以测站坐标作为未知数。由于采用双差观测值作为基本观测量，在测站之间按全组合形成不同的基线。GAMIT是将所有基线的观测方程一起处理，只建立一个法方程，一次性解算出所有未知点的坐标，在O文件中以基线形式输出。无论精度如何，同步环闭合差总为0。在O文件中，基线形式以直角坐标系和站坐标系两种形式给出，即(DX，DY，DZ，S)和(N，E，U，S)以及各个分量的标准差。站坐标与直角坐标的转换公式为：

需要说明的是，基线分量的协方差矩阵是根据站坐标未知数的协方差(法方程系数阵的逆和观测值单位权方差相乘)计算的，NEU分量的标准差可以由(DX，DY，DZ)的方差-协方差矩阵根据误差传播律计算。在基线较短时，一般将站坐标系中基线NEU分量的误差作为基线水平方向和大地高方向的误差；基线较长时应考虑基线NEU分量的精度与测站点NEU分量的精度之间的差别。如图3-2所示，设基线方向为测站A到测站B，站坐标系的原点为A点，假设A站点坐标没有误差，基线分量的误差即代表B站点的坐标误差。基线较长时，U方向的误差几乎是B点水平方向的误差。因此，在分析测站水平和大地高方向的误差时，如果基线较长，不能仅从NEU基线结果的精度来分析，而应以B站的站坐标精度为准。

图3-2　测站误差分析示意图

以拉萨至上海基线为例，基线长为2865km，当基线NEU方向变化(0，0，40)mm时，反映到上海站站坐标的位移量是(1，-17，36)mm。在分析测站水平和大地高方向位移或精度时，如果以基线为对象应考虑这一微小区别。

(3)数据处理质量的评价标准

GAMIT基线解算结果的好坏，一般有以下几种评判标准：

①GAMIT解算结果中的标准化均方根误差NRMS(Normalized Root Mean Square)用来表示单时段解算出的基线值偏离其加权平均值的程度，是从历元的模糊度解算中得出的残差。NRMS是衡量GAMIT解算结果的一个重要指标，其计算公式如下：

(www.xing528.com)

式中，Yi为基线向量历元解算值，Y为基线向量真值，N为历元总数，为各历元解算值中误差。

一般说来，NRMS值越小，基线估算精度越高；反之，精度较低。其值一般应小于0.3，若NRMS太大，则说明处理过程中周跳可能未得到完全修复，或存在其他问题。

②参数的改正量不能大于其约束量的2倍。

③当Choice of observable为L1_ONLY时，B1Ll计算的整周模糊度必须是整数。

④一般以坐标的重复性作为衡量坐标解算结果的评价指标。Xij，Yij，Zij表示j点在i测段(i＝1，2，…，n为测段数)算得的坐标，则坐标分量重复性为：

式中，σXj，σYj，σZj分别为点的坐标分量重复性，PXi，PYi，PZi为i测段解的坐标分量的中误差平方倒数，为坐标分量加权平均值，可分别由式(3-22)求得：

⑤基线重复性是衡量数据处理质量的重要指标之一(刘经南，等，1995)，可采用以下两式分别计算基线向量的重复性和相对重复性：

式中，Rl为基线向量的重复性，Rr为基线向量的相对重复性，n为基线单日解数目，Li为第i日的基线分量(或边长)，为单天解基线分量(或边长)的加权平均值，其公式如下：

进一步以基线重复性为观测值，用线性拟合求出重复性的常数部分和与边长成比例的部分：

2.GLOBK软件功能及组成

(1)GLOBK平差方法

GLOBK平差方法是基于动态参数的估算方法——卡尔曼滤波，要求观测方程中某些参数满足某一动态方程。卡尔曼滤波的数学模型包括观测方程和状态方程，如式(3-27)，其中，εk表示观测噪声，ωk-1表示动态噪声，Lk为观测值，Xk为状态值，Ak，φk分别为观测方程和状态方程的系数阵。

卡尔曼滤波的随机模型如下式所示：

式中，Pε、Qω分别表示动态噪声和观测噪声的方差矩阵；δ(tj-tk)表示kronecker函数，该函数满足的条件为：

系统的初始状态X0是一个随机向量，该向量具有正态分布或其他分布，其均值和方差阵为：

满足式(3-28)，式(3-30)的卡尔曼滤波递推公式可表示为：

对于两期的观测结果，经过相似变换可获得两期的ITRF坐标，两期坐标作差便可以得到ITRF框架下测站的位移，除以时间则为测站速度；对于多期观测，一般将参考时刻的坐标和速度当作未知数建立方程求解，对GAMIT基线解算所获得的准观测值，建立如下观测方程：

式中，ΔX0表示历元t时刻的松弛解坐标；X0、ΔV0分别代表参考历元t0的坐标和速度；λ为尺度因子；τX和τV为参考历元的平移参数和变化速率参数；ωX和ωV为参考历元旋转参数和变化速率参数；γ(t)为随机噪声；δεk为第tk时刻的位移突变，在估算测站速率的时候应予以考虑，并有

且μ有如下的表示形式：

参考时刻的坐标先验值用X0(x0，y0，z0)表示，待估参数用p表示，即p＝(X0，ΔV，λ，δεk，τX，τV，ωX，ωV)，那么状态方程可以表示为：

根据观测方程式(3-36)和状态方程式(3-39)，可按经典的卡尔曼滤波估计参数。

(2)GLOBK软件模块

GLOBK软件模块可大致划为四大类：

1)格式转换模块(htoglb)

这个模块是将由GPS、VLBI和SLR等分析软件的解文件转换成GLOBK软件所需要的二进制文件h-文件。目前支持如下几类文件：

①GAMIT软件的h-文件；

②GPS(或其他空间大地测量技术)SINEX格式文件；

③FONDA软件的h-文件；

④JPL机构提供的Stacov文件；

⑤包含站坐标和速度场的SLR/GSFC文件；

⑥包含站坐标和速度场的VLBI/GSFC文件。

2)运算模块(GLRED，GLOBK，GLORG)

GLRED模块通过调用GLOBK模块分析单天解，其生成的解文件可以用来形成时间序列。GLOBK模块是GLOBK软件的主模块，实现该软件的功能。GLORG模块可以为平差结果定义参考框架，具体通过固定(或约束)站坐标和速度由坐标转换来实现。GLORG模块可以单独运行，也可以被GLOBK/GLRED模块调用。所涉及的重要文件是cmd_file，它们是GLOBK和GLORG的控制文件，内容包含求解策略等。

3)GMT图形应用模块

这类模块主要包括：sh_plotcrd、sh_globk_scatter、multibase、sh_plotvel等，主要功能是利用GMT软件绘制时间序列、速度场等图形，可用于分析数据质量和测站的地壳运动等情况。

4)其他辅助模块

其他辅助模块主要包括：glist、glsave、extract、exbrk、corcom、cvframe、velrot等。这里面有两类，一类是为GLOBK等模块服务的，如glist、glsave等；一类是用于框架之间和板块运动分析的，如corcom、cvframe、velrot等。

(3)GLOBK常用文件及格式说明

1)输入文件：GPS、SLR、VLBI和SINEX文件

随机特征可由apr_XXX和mar_XXX表述。

2)数据文件和控制文件

包括二进制H-文件和指令(cmd)文件。

3)GLOBK结果文件中的NEU坐标

GLOBK的输出结果文件一般为∗.prt和∗.org，在给出ITRF坐标的同时还给出了新的NEU坐标，它与站坐标定义的NEU不同，这种坐标类似于平面坐标，属于圆锥投影。由X、Y、Z先计算出测点的大地经纬度和大地高U，直角坐标与大地坐标的转换公式为：

式中，ae为WGS84椭球的长半径，e2为第一偏心率。N、E应严格定义为：

式中，B、L的单位为弧度，r0为余纬为θ0时的纬圈半径，即

θ0定义为最接近0.00005弧度(约10角秒)的余纬，

如此定义的NEU坐标意义在于：经度方向的平差值不受纬度方向的微小变化而变化，而NEU方向的中误差仍以站坐标的形式表示。

4)运行GLOBK时的注意事项

①GLOBK是基于线性模型的。因此，在测站坐标或轨道参数的改正值较大时(测站坐标改正值大于10m或轨道参数的改正值大于100m)，需要进行前期数据的再处理以获得满足要求的准观测数据。

②GLOBK不能解决在前期数据处理阶段因周跳未得到完全探测、数据质量差或大气层延迟模型误差所带来的问题。在GLOBK数据处理阶段，不能彻底消除特定测站或卫星的影响，只能通过特定手段减弱其影响。

③GLOBK不能进行整周模糊度的解算，因此，在前期数据处理阶段必须完成整周模糊度的解算。