一种航空重力测量GPS后处理方法与流程

本发明专利涉及勘探地球物理学，GPS资料在航空重力探测中应用，即利用GPS资料高精度后处理解算飞机的位置、速度和加速度等。

背景技术：

地球重力场测量和研究对大地测量、军事科学、空间技术、海洋学等学科都发挥着举足轻重的作用；此外，它也是检测火山喷发、局部地层沉降、抗震减灾的重要手段。航空重力测量是以飞机等为载体的一种新型的重力测量技术，除具有快速、经济、灵活、精度高等特点外，而且能够在一些难以进行地面重力测量的区域进行作业，如冰川、沼泽、原始森林、沙漠等交通不便和无人居住区。

GPS在航空重力测量中不可或缺，这是因为确定了载体飞机的位置资料，重力测量才有意义；有了速度、高度等参数才能厄特弗斯改正、上下延拓处理等；另外，载体的加速度是航空重力测量值中的噪声，需要从观测值中去除，而这些问题的解决需要利用GPS技术来解决。航空重力测量系统本身也是集多个分系统、多种技术于一体的新型重力测量装备，而放置飞机内部的重力传感器（比力仪）测量数据中无法区分属于飞机运动而产生的扰动加速度还是地下不均匀地质体产生的重力加速度，因此关键技术之一是从重力仪所测总加速度中分离出飞机运动等因素产生的扰动加速度。精确的确定扰动加速度依赖于精确获得载体的位置、速度和加速度等飞行状态参数，并利用这些参数对重力比力仪所测结果进行必须的改正。

目前，GPS的实时处理、后处理方法和软件较多，但主要解算GPS定位（位置解算）的问题。因此，需要GPS后处理方法来获得航空重力测量中飞机的位置、速度和加速度的方法。

技术实现要素：

本发明的目的在于：提供一种基于GPS原始测量资料，采用后处理方法，高精度解算航空重力测量的载体的时间、位置、速度和加速度的技术方法或方案，数学表示为或，称之载体十维动态向量（或十要素），而解决现有的GPS处理方法无法直接满足航空重力测量的技术要求。

为解决上述技术问题，本发明采用以下技术方案。

航空重力测量系统中的GPS子系统包括机载GPS动态接收机和地面GPS基站接收机两部分，距离一般小于5km，主要用于获取GPS授时、GPS卫星广播星历、GPS观测数据等信息，可利用GPS接收机的授时信号进行不同传感器之间的时间同步。其中地面GPS接收机为基站，即固定位置静态测量；机载GPS接收机为动态测量，测量频率为1Hz～10Hz，即1s～0.1s测量一次数据。此外，飞机起飞前需要静态测量不少于15分钟，为获取GPS相位测量的整周模糊度，同时也为航空重力测量取得地面的重力基准值。测量的原始数据需转换为与接收机无关的标准数据交换格式——RINEX格式。

采用后处理方法步骤如下：

（1）读取基站观测文件、机载GPS观测文件、广播星历文件或下载精密星历文件；

（2）利用星历文件或精密星历文件计算GPS卫星在WGS-84坐标系下的历元、位置、速度等状态参数；

（3）由于卫星广播星历文件每小时更新一次，精密星历文件15分钟更新一次，而GPS接收机0.1s或1s观测一次，因此将上述的卫星结算结果采用插值方法计算与观测文件历元序列（时间序列）一致的卫星运行的轨道参数，按卫星编号输出与观测数据的历元序列（时间序列）一致的独立文件或数据表；

（4）地面的GPS基站位置为已知数，计算不同历元与观测得GPS卫星之间的空间几何关系，即卫星与基站之间的距离、方位角和高度角，因为高度角小于20°的卫星计算误差较大，故舍弃观测的高度角小于20°的卫星；

（5）若观测的卫星数目不少于4颗时，可利用机载GPS动态接收机观测文件的伪距、相位、多普勒频移等观测量计算飞机时间、位置、速度和加速度（10维状态向量），采用相位计算飞机位置时，需要计算整周模糊度。其计算方法可采用最小二乘法，即GPS单点计算方法；或者

（6）对地面GPS基站和机载GPS动态接收机的同一历元观测的同一卫星的伪距、相位和多普勒频移做差分处理，即单差分处理；在单差分基础上，观测的不同卫星进行再次差分处理，即双差分处理；同样采用最小二乘法计算GPS动态接收机的时间、位置、速度和加速度（10维状态向量）；

（7）由于GPS接收机天线的相位中心和重力传感器重心并不一致，需将GPS天线的位置归一到重力测量仪中心位置，即偏心改正处理；

（8）应用数值计算方法高精度获得飞机“真实”飞行姿态参数，方法包括滤波器、最优估计方法的应用、数值信息处理方法等。

本发明的有益效果在于：通过对GPS测量的原始资料后处理解算，获得飞机飞行十维状态向量参数，即GPS时间、位置、速度和加速度等飞行姿态参数，特别是利用GPS技术使得航空重力测量总加速度中分离出飞机运动等因素产生的扰动加速度。此外，GPS后处理技术可广泛的应用于空间大地测量、航空摄影测量和全球性的板块运动、地壳运动检测、地质灾害监测等领域。

附图说明

下面结合附图和实施例对本发明进一步说明。

图1，航空重力测量中GPS十维动态向量定位解算流程。

图2，导航文件计算GPS卫星轨道参数及数据标准化流程。

图3，观测文件数据标准化流程。

图4，计算位置、速度、加速度的流程。

图5，导航文件计算GPS卫星状态实例。

图6，解算的飞机的位置、和速度实例。

图7，解算东向加速度曲线实例。

图8，解算北向加速度曲线实例。

图9，解算垂向加速度曲线实例。

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行进一步说明和描述。显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。

如图1所示，航空重力测量中载体十维动态向量解算流程原始数据包括GPS的观测数据文件和导航星历文件，其中导航文件用于提供解算GPS卫星轨道的资料，观测数据文件用于提供解算GPS接收机的时间、位置、速度和加速度（即本发明所述的十维动态向量）。若仅包括移动站观测数据文件和导航文件，则为单点十维动态向量解算；若包括移动站和地面基站观测数据，则可高精度的差分十维动态向量解算。

如图1所示，读取了导航数据后，首先需要解算观测历元时刻的卫星位置、速度等。一般而言，导航文件每四小时更新一次，而观测数据的采样间隔为0.1秒或1秒，在此期间的卫星相对于地球不间断的动态运动中，而解算载体的位置和速度需要历元相对应的卫星轨道参数，换言之，接收机观测数据和导航文件提供的卫星轨道参数历元时刻可能不一一对应。因此，根据卫星的运动轨道理论，将利用星历文件确定的卫星轨道参数进行内插获得任意历元（一般为观测历元相应）的卫星位置和速度，称为轨道标准化。

卫星的轨道参数常作为时间多项式，根据已知的卫星星历文件提供的信息计算指定时刻的卫星位置和速度；将观测历元变换成，即，；则卫星坐标、速度和加速度的分量可用如下多项式表示：

，，

其中和分别为开始历元和拟合时间区间的长度，为多项式的阶数，为多项式系数，它的导数用递推方式确定：

然后根据已知的卫星轨道状态参数，利用最小二乘法可以求出；进而确定不同历元时刻的卫星状态参数。

计算观测的某一个卫星的高度角、方位角和距离的随历元变化，其解算流程如图2所示。

如图1所示，将卫星轨道的历元序列数据格式和观测数据的历元序列数据格式进行历元同步和数据融合，进而计算接收机与观测卫星之间的几何关系，由于卫星观测角度对计算误差的影响不可忽视，故舍弃观测高度角度小于20°的卫星。

如图1所示，要求移动站观测不少于4颗卫星的数据，利用观测的伪距、相位参量，采用最小二乘法计算移动站的位置，利用观测的多普勒频移动计算移动站的速度，即所谓的GPS单点解算方法。将基站和移动站观测的同一个卫星的数据做差分处理，可消除该卫星的钟差影响和部分误差，即单差分计算方法；在单差分基础上，对观测不同卫星的数据再次差分，可进一步消除传播路径和接收机部分误差，即所谓的双差分计算方法。解算流程如图3所示。

实际的航空重力测量中，为不影响GPS卫星信号的接收，GPS天线一般安装在飞机顶部，而重力仪安装在机舱的底部，显然重力传感器中心与GPS天线相位中心并不一致，而航空重力各参量计算都应归化到重力比力仪的摆杆中心，因而GPS天线相位中心亦应该归一至重力比力仪的中心。若飞机平稳飞行时，GPS测量位置加上相应的距离即可改正为重力仪的位置，GPS测量的速度、方位和加速度即为重力仪的速度、方位和加速度；然而，由于飞机飞行的俯仰和横滚，GPS测算结果归算过程中必须考虑飞机姿态对GPS解算的影响，这种由于GPS中心位置偏离重力仪中心的位置的影响的必须改正。在归算过程中，因为GPS测量结果在地固坐标系，测定的姿态参数、重力测量信号均在不同坐标系，因此，需要把GPS测量结果化至相应坐标系。改正可以应用如下矩阵变换方法： GPS数据的测量结果属于地固坐标系(WGS-84)，首先需要将其转换至水平坐标系，位置转换根据如下转换公式：

其中：

式中， 、、分别为飞机的航高、经度和纬度， 为地球卯酉圈的曲率半径，为地球椭球的长半轴，为椭球第一偏心率。对上述公式式两边求一阶导数，可得速度坐标变换公式；同理，对上述公式求二阶导数可得加速度换算公式。 航空重力测量中，以飞机为载体，飞机重心为原点，轴为通过的机身中轴线，指向机头方向为正；轴为过的机翼中轴线，指向右为正；轴过垂直向下为正。建立载体直角坐标系，设、分别为重力传感器中心和GPS天线相位中心在中的位置矢量；为为重力传感器中心和GPS天线相位中心在中的位置矢量差，将转换至当地水平坐标系公式如下：

其中：

式中，、、分别为载体的横滚角、俯仰角和方位角，一般统称为姿态角。改正解算的GPS天线的相位中心归一重力仪的重心，即可得到重力仪器重心的不同历元的位置、速度和加速度等状态参数。

根据如图1所示，对速度矢量序列进行差商计算，即可获得加速度矢量序列。差商计算方法如下：

从而获得飞机对重力传感器产生的扰动加速度 。

采用上述流程，可以利用GPS资料获得航空重力测量中的飞机时间、位置、速度和加速度十维状态。

但由于误差的和干扰因素的存在，解算结果往往含有噪声，需采用数值方法消除或减弱误差影响。

如图1所示的滤波器包括时间滤波器、空间滤波器和二维滤波器。首先，航空重力测量GPS测量是与时间密切有关，飞行的噪声频率往往较高，而重力有效信号属于低频信号，故采用时间低通滤波，从而消除或减弱高频干扰。理想低通滤波器是一个过渡带为“零”的精确矩形函数，其通带的函数为“1”，阻带的函数值为“0”，这在物理上或者数字信号处理中是无法实现的。因此，可实现的低通滤波器是理想低通滤波器的一种逼近。一个长度为的FIR滤波器的频率响应可表示为：

根据数字信号处理的原理，可写为模和幅角的形式：

上式中是的增益()，是一个可正可负的实函数，是相频特性。为了便于实现，可将化为一个ω的固定函数和一个余弦求和式的乘积，即：

假设被逼近的理想频率响应为，则逼近的加权误差为：

其中是权函数，表示通带和阻带具有不同的逼近精度：

于是，用一致逼近的问题可描述为：求系数组使在被逼近的频带范围内误差函数的极大绝对值为极小。由交错点组定理及Remz算法可解决该逼近问题，从而完成滤波的设计。

球物理测量中，通常把物理场的信号看做空间位置的函数。因此，也可依据地球物理场的特征将的GPS动态测量信息看做空间位置的函数。因此，在航空重力测量中的GPS解算中，也可把GPS解算结果（垂直加速度）作为空间位置的函数。设计空间域滤波器和时间域滤波器一样，都是一维滤波器的设计；故涉及的理论和数学方法也相同。地球物理场的信号处理中的位置大都是三维空间的位置关系。在设计滤波器中，可利用空间位置的某一点作为基准点，利用空间位置与该点的距离，然后按信号与的函数关系：

上式中，

可以利用同样的滤波器。根据一个基准点，由于空间位置可能有三个点与距离相等，因此，在实际中可采用两种处理方法：第一种就是按类似GPS定位的方法，可设计3个或3个以上的基准点，然后的定位方法解算空间唯一的位置；另外一种方法就是分段（有唯一的，如按直线）解算。一般来说，在空间域，重力场异常的有效波频率在0.01Hz以下（或周期100s以上），而飞机的垂向加速度（干扰波）则在0.1Hz以上（或周期10S以下）；故利用空间滤波器基本能够分离出重力异常的有效信号。

在实际的航空重力测量中，无论是GPS解算出加速度，或是重力（比力）仪器测量的物理场，可能有效波和干扰波的频谱成分十分接近，甚至重合，这时无法利用空间域或时间域频率滤波压制干扰。因此，需要利用有效波和干扰波在其他方面的差异来进行滤波，这种滤波器就是二维滤波器的应用。而航空重力的测量实际是时间和空间的二维函数，因此，即可用振动图来描述，也可用波剖面来说明它，二者之间连续如下公式：

上式中为空间波数，表示单位长度上波长的个数；为频率，描述单位时间内振动次数；为波速。显然航空重力测量的二维函数与飞机飞行的速度有关，并且航空重力测量总是沿测线进行观测，相对某个基点，上述的波数和速度应以波数分量和视速度代入；则由如下公式：

既然航空重力信号（加速度）不是单独的空间的函数，也不是单独的时间的函数，且空间和时间存在着密切的关系，无论单独进行哪一维的滤波都会引起另一维特性的变化，例如单独进行频率滤波会改变剖面形状，单独进行波速滤波会影响振动图形，产生频率畸变，即单独采用时间滤波或空间滤波可能会产生信号特性畸变，产生不良效果，那么只有根据二者的内在联系组成时间空间域（或频率波速域）滤波，才能有效压制干扰、突出有效波的目的。二维滤波是建立在二维傅里叶变换基础上的。而航空重力测量的加速度是一个随时间和空间变换的波，通过二维正、反傅里叶变换得到其频率波速谱和频波谱的时空函数，如下公式描述：

上式说明，是一个无数个频率为、波数为的平面简谐波所组成，它们沿测线以视速度传播。二维滤波的计算可根据二维滤波原理（二维傅里叶变换）进一步推导。从而把时间-空间域的函数关系转换为频率-波数域的函数。如果有效波和干扰波的平面简谐波成分有差异，有效波的成分以与干扰波的平面简谐波成分不同的视速度传播，则可用二维视速度滤波将它们分开，达到压制干扰、提高信噪比的目的。为了有效利用信号的二维函数特性，有两种方法，即：同一测线，利用不同的飞机飞行速度反复测量（即不同），获得不同时空函数关系，从而获得信号的二维特征。另外就是把测区作为整体，选择合适的滤波参数，利用二维函数。同样道理，航空重力解算的加速度数据可看作时间和位置相关函数，故可采用二维滤波器消除或减弱干扰因素的影响。

如图1所示的流程图，在航空重力处理解算中，为了更好的逼近飞机飞行的真实状态，采用最小二乘法将解算的结果最优估计处理，建立描述飞机10维状态的拟合数学模型。最小二乘法其基本思想是残差的平方和为最小，对数据，求次多项式：

使总误差的平方最小，即：

把看作关于的多元函数，故上述拟合多项式的构造问题可归结为多元函数的极值问题，并以此建立方程组：

利用最小乘法可以最优估计飞机的10维状态参数。

如图1所示的流程图，在航空重力处理解算中，为了更好的逼近飞机飞行的真实状态，也可采用Kalman最优估计的方法将解算的结果最优估计处理。在GPS动态定位中，运动接收机天线遵循一定规律不断变化的动态系统，不同时刻的系统状态与所进行的观测存在着相互联系。对状态估计时，这些规律和联系应当充分加以利用。对于离散时间系统，状态随时间变化的规律通常用向量微分方程来描述，称为动态方程或状态方程。这个方程受到某种随机扰动的影响，这种扰动称为动态噪声，观测向量与状态向量或动态向量之间，存在着一定的函数关系，称之为观测方程；同时，还存在着随机观测误差，称为观测噪声。Kalman滤波正是充分的利用了上述规律的最优估计理论的计算方法。利用kalman滤波器处理动态数据，首先需要建立滤波的动态模型（状态方程）和观测模型（观测方程）。常速模型和常加速模型的选择依赖于运动载体的运动状态以及数据采样率的高低，在高精度GPS动态定位中数据采样周期一般为1s或更小，此时可采用常速模型。动态噪声可假定为零均值的高斯噪声，选取状态向量，其中分别为接收天线在WGS-84地心坐标系中的三维位置，以弧度为单位，以米为单位；，，分别为动态GPS接收天线北向、东向、垂向的速度；，分别为GPS接收机的时钟偏移和时钟漂移。则利用Kalman滤波处理GPS动态向量定位数据的滤波模型的状态方程为：

观测方程：

上式中，为观测值；为转换矩阵；分别为、时刻的量测值；为干扰矩阵；为状态噪声；为观测噪声；为观测矩阵。即Kalman滤波器的输入量为GPS解算结果，取观测方程如下：

为固定零均值干扰的高斯任意序列，采用高斯任意序列方程，在飞机飞行时，由于速度与位置相关，其自相关方程如下：

其中：是自相关值；是相关时间的倒数(1/s)；是系统噪声的标准偏差。为了使用Kalman状态方程，我们首先要令未知参数和作为参数的函数。对上式两边同乘以时刻的状态函数，得如下方程：

假定不相关且：

所以，

故可得如下方程：

对状态变量进行平方，可得如下的期望值：

由于,参数，根据,，那么可得如下方程：

其中包括以下状态变量：纬度（latitude[rads]）；经度（longitude[rads]）；高度（height）[rads];北向速度m/s;东向速度m/s;向上速度m/s;接收机时钟偏差m；接收机时钟漂移m/s。

为了进一步验证本方案所达到的效果，下面给出本方案在实际应用中的实施例 ：

实际数据来自某次航空重力测量的数据，测量的起止时间为2007年10月02日23:00:16.000和2007年10月03日03:05:29.000，采样间隔为1s，飞机起飞前后有10～15 min的静校正时间，利用L1相位载波进行解算。

图5是根据该航空测量的GPS资料的星历文件计算了卫星位置、速度，且经过插值结果。其星历文件更新频率为1个小时，数据标准化后的卫星状态参数为1s，与观测文件的采样频率一致。

图6输出的解算位置、速度结果，经滤波、最优估计后的结果，曲线相对圆滑、噪音消弱，更符合飞机真实的状况。

图7~图9分别取上述解算的速度结果，经偏心改正后，采用差商方法解算的飞机东向、北向和垂直加速度状态曲线，可清晰分析飞机飞行中对重力加速度传感器起在不同方位的干扰。 尽管为了更逼近的获得“真实值”，可使用各种处理手段，由于个人经验和偏好的不同，选择参数和处理流程也存在差异，也可能导致解算结果略有出入。需要根据实际资料分析，甚至反复试验而获得合适的处理参数。

以上仅是本发明的优选实施方式，本发明的保护范围并不仅局限于上述实施例，凡属于本发明思路下的技术方案均属于本发明的保护范围。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理前提下的若干改进和润饰，应视为本发明的保护范围。