1 引言

我国重力空白区域在我国国土面积中占有很大比重，并且多集中在难以到达的沙漠、水域、高山等区域(陈俊勇等，2001；周坚鑫等，2003；张昌达，2005).传统的静态重力测量手段难以在这些困难区域展开，卫星重力测量由于分辨率较低，也难以满足需求.而航空重力测量以其高精度、高分辨率、高机动性等特点，成为该类地区获取重力信息的有效手段.航空重力测量的研究始于20世纪50—60年代，在80年代末90年代初逐渐成熟.研究表明航空重力测量的精度为2~10 mGal(分辨率6~10 km)(Brozena et al.，1989；Klingelé et al.，1995；Brozena and Childers，2000；Forsberg et al.，2001； Verdun et al.，2002；Rene et al.，2007；孙中苗等，2013a).商用航空重力仪(例如LaCoste & Romberg、BGM、KSS等产品)具有较高的精度和较好的稳定性，但我国对该类设备的引进存在诸多限制.我国在航空重力测量的研究起步较晚，西安测绘研究所和航天科技集团16所率先研制的航空重力测量系统CHAGS，测量精度达到了2~5 mGal(分辨率10 km)(夏哲仁，2004；孙中苗，2004；孙中苗等，2004a，2004b；孙中苗等，2013a).

基于GNSS/INS的航空重力测量是目前研究最为广泛的航空重力测量技术，根据测量原理的不同，可以分为直接求差法和状态模型法(Jekeli，2001).直接求差法由加拿大Calgary大学的Bruton在2000年提出(Bruton，2000)，随后Serpas、Senobar等对该方法进行了深入的研究(Serpas and Jekeli，2005；Senobari，2010).在国内，孙中苗等利用直接求差法对CHAGS系统开展了多项研究(孙中苗，2004；孙中苗等，2013b).尽管直接求差法在航空重力测量领域中有较为深入的研究，但该方法对惯性系统的校准精度和系统姿态的测量精度都有严格的要求.与直接求差法不同，状态模型法把重力扰动当作Kalman滤波系统的状态量进行最优估计，在一定程度上降低了对惯性系统的要求.Tomé、Deurloo等的研究成果表明，利用中低精度的惯性系统，根据状态模型法同样可以得到较高精度的重力扰动结果(Tomé，2002；Deurloo et al.，2011；Bastos et al.，2013).目前国内的研究多集中在直接求差法，对状态模型法的研究较少.

自适应Kalman滤波以其调节状态信息的权阵来控制状态误差对参数估计的影响的特性，在定位与导航有着广泛的研究与应用(Yang et al.，2001，2006a，2006b；Hajiyev and Soken，2013).研究表明引入自适应因子，并根据一些理论或者经验的准则对自适应因子进行合理的估计以后，组合导航系统的性能明显优于经典Kalman滤波(高广为等，2006；Yang et al.，2006a).

本文首先展开对状态模型法航空重力理论研究，给出该方法详细的实现过程，并在状态模型法的基础上引入最优自适应因子算法，调节动力学模型信息和观测信息对滤波结果的影响，改善重力扰动的计算结果，并以新疆地区的实测数据验证该方法的有效性.

2 基本原理

根据牛顿第二定律，在惯性坐标系中质点的动力学方程为

其中 为质点的惯性加速度，gⁱ为重力加速度，fⁱ为比力加速度.

在SINS/GNSS组合导航过程中，在导航坐标系n下，速度测量值对时间求偏导可以得到(Titterton and Weston，2004)：

其中v_eⁿ为载体的速度，f^b为加速度计的比力测量值，C_n^b为b系到n系的方向余弦矩阵，ω_ieⁿ为地球自转角速度在n系下的投影，ω_enⁿ为n系相对于地球坐标系e系的角速度在n系下的投影，式中的第二项为哥式加速度和向心加速度在n下的投影，gⁿ为实际重力矢量.gⁿ通常采用正常重力γⁿ和重力扰动δgⁿ的矢量和表示，公式为

在物理大地测量的定义中，重力异常主要包括重力扰动(纯重力异常)和混合重力异常.其中重力扰动是指地球某一点上的重力观测值和正常重力之差.重力扰动δgⁿ为一个三维矢量，由两个水平分量ζg，ηg和一个垂直向下分量δg组成(Heiskanen and Moritz，1967)，δgⁿ=[－ζg －ηg δg]^T.本文使用的惯性单元为中等精度水平的光纤陀螺惯导，受惯性单元性能的影响，δgⁿ水平分量很难达到期望的精度，因此这里主要讨论垂直分量上的重力扰动δg.目前估计重力扰动δg的方法有两种：直接求差法和基于Kalman滤波的状态估计法，下面对两种方法的基本原理进行介绍.

2.1 直接求差法航空重力测量

直接求差法主要是通过公式将重力扰动矢量表示为

重力扰动矢量垂直向下的分量可以详细的表示为

其中(·)_D为垂直向下的分量，右边第三项为厄特弗斯改正，L为纬度，h为椭球高，v_E为东向速度，为垂直向下加速度，R_M为子午圈半径，R_N为卯酉圈半径.

利用公式(2)计算的空中重力扰动含有大量的噪声，而重力扰动信号能量通常集中在低频段很窄的部分.因此需要设计合适的低通滤波器对空中重力扰动进行噪声消除.在航空重力测量领域中，比较典型的低通滤波器有RC低通滤波器、Butterworth低通滤波器、FIR低通滤波器.FIR低通滤波器的差分方程形式可以表示为(孙中苗，2004；孙中苗等，2004b)

其中x(n)为输入信号，y(n)为输出信号，h(k)为滤波器系数(单位冲激响应)，N为滤波器长度.h(k)的理想频率响应可以表示为

确定单位冲激响应函数主要有切比学夫等纹逼近法和窗函数法，主要做法请参考文献(孙中苗，2004)和(孙中苗等，2004b)，本文不再赘述.

2.2 引入最优自适应因子的状态模型法航空重力测量

在组合导航领域，位置、速度、姿态、惯性器件的系统误差通常被当作Kalman滤波的系统状态量，通过预计方程和更新方程估计出来.基于状态法的航空重力测量就是通过对重力扰动进行建模，将其作为系统的状态量，利用Kalman滤波进行最优化估计.

以21+1维系统状态的Kalman滤波为例，系统的状态量可以表示为

其中r_eⁿ、v_eⁿ分别为位置向量和速度向量，C_bⁿ为从载体坐标系b旋转到导航坐标系n的方向余弦矩阵，b_a和b_ω分别为惯性单元加速度计和陀螺三轴方向的随机常值误差，s_a和s_ω分别为惯性单元加速度计和陀螺三轴方向的尺度因子.δg为重力扰动在垂直向下方向上的分量，假设在短时间内为一常量偏差，也可以采用三阶Gauss-Markov模型进行建模(Jekeli，2001)，研究表明两种模型都能很好的对重力扰动进行估计(Deurloo，2011).对应的误差状态向量可以表示为

Kalman滤波的动力学预计方程为： =f(x)+Gw，其中w为过程噪声，主要由加速度计和陀螺观测量的随机噪声组成.利用一阶泰勒公式进一步展开，可以得到误差状态量的预计方程的近似表达方式： ，其中F是系统矩阵，其非线性的表达方式为 ，G和G′为噪声影响矩阵，F和G都是时变矩阵，在本文的应用中表示为

其中F₁₁—F₃₃由位置、速度、姿态的导航方程求偏导得出，读者可以参考相关文献或自行推导；C_bⁿ为从惯性系到导航系的方向余弦矩阵.

一个典型的Kalman滤波的预测方程可以表示为

其中Φ_k－1近似的公式为

噪声矩阵Q_k－1为：Q_k－1≈G_k－1QΔtG_k－1^T.其中功率谱密度矩阵Q可以由过程噪声表示为(Gelb，1974): Q(t)δ(t－τ)=E[w(t)w(t)^T]，在本文的22维的Kalman滤波中，Q被定义为

其中Q_a，Q_ω分别是加速度计和陀螺随机游走的平方，q_δg=σ_δg²，这里σ_δg取 (Bruton，2000).

引入最优自适应因子的更新方程可以表示为

其中α为最优自适应因子，z_k是第k时刻的更新信息，通常是指GNSS的位置和速度的解算结果，H_k为观测矩阵，Kalman增益矩阵K_k为

其中R_k为更新观测值的协方差矩阵.

t_k 时刻预测状态向量的协方差可以表示为

预测向量的不符值为： ，利用Sage开窗估计法，可以求得状态预测向量协方差矩阵的新的估值为(杨元喜，2006)

采用多个历元的平均值并不能有效的反应当前历元的模型误差的大小，因此这里简单的取：

自适应因子的表达式为(Yang and He，2001；高为广等，2006；Yang et al.，2006a，2006b；杨元喜，2006)

其中 ，c在本文中取1.

3 算例及结果分析

选取西部某地两天航空摄影测量任务的POS(测姿定位系统，Position and Orientation System)数据进行分析，飞行轨迹分别为图 1和图 2.搭载的POS系统为IGI AeroControl & AeroOffice，该系统包括一款中等精度的光纤陀螺惯导IMU-d和一款NovAtel OEM4 GNSS双频接收机，表 1为该系统利用本文的算法计算的位置和姿态结果与航空摄影测量的空中三角测量结果的较差统计表，从表可以看出利用该系统能提供比较可靠的位置和姿态信息.飞行载体的为运八型飞机.在测区内架设有GPS基站，利用GPS差分技术获取高精度的位置、姿态、速度、加速度信息.测试过程中IMU设备的频率为256 Hz，GNSS差分结果的频率为1 Hz.分别利用直接求差法和自适应状态模型法进行航空重力解算.解算结果与全球重力场模型EGM2008在该区域重力扰动资料进行比较，验证结果的精度.EGM2008全球重力场模型是由NGA(National Geospatial-intelligence Agency)发放全球超高阶地球重力场模型，分辨率最高可达2.5′×2.5′(约为5 km格网)，研究表明该模型在我国地势平缓的区域的精度能达到10 mGal左右(章传银等，2009；杨金玉等，2012).由于目前掌握该测区的地面重力资料非常有限，本文暂时用EGM2008全球重力模型计算出的重力扰动值作为参考，对初步解算结果进行对比分析.两次测量都在飞机起飞前进行了系统初始对准，系统的初始对准包括POS系统的初始姿态对准和重力扰动的对准，对准时间约为15 min.由于停机坪附近没有重力基准点，系统采用EGM2008模型推算的重力扰动作为初始重力扰动.EGM2008提供的全球重力扰动结果为参考椭球表面的重力扰动，本文利用NGA网站提供的Fortran脚本“EGMhsynth_WGS84.f”和2.5′×2.5′的格网数据进行计算了测试所在区域椭球表面的重力扰动，然后再通过向上延拓公式延拓到航线所在的高度上，将延拓后的结果作为本文计算结果的参考值.测试中采用的向上延拓算法为直接代表法，计算公式为(石磐和王兴涛，1997；孙中苗，2004)

其中δg_EGM为EGM2008计算的椭球表面的重力扰动，b= 0.1119 mGal·m^-1为布格系数，H_T为航线高度，为平均航线高度.

在k时刻的重力输入值的迭代公式为：g_k^air=g_k^h+δg_k-1^air，然后利用公式(8)—(15)计算ε_δg^k，k时刻的重力扰动就可以表示为：δg_k^air=δg_k-1^air+ε_δg^k.其中g^h为在航线高度上的正常重力，公式为(NIMA，2000)：

其中

这里的参数Φ为纬度，单位为弧度，h为椭球高，单位为m；a，b，e，f是重力模型的椭球参数；参数γ_p和γ_e分别为重力模型的极参数和赤道参数；ω_e为地球的自转角速度，GM为地球的引力常数.这些地球参数的取值请参照文献(NIMA，2000).

3.1 算例1

算例1为对新疆某沙漠地带一平均航高6000 m、航迹340 km、航速155 m·s^-1的航线进行航空重力解算，飞行轨迹如图 3.飞机在飞行过程中受到大气扰动较小，飞行高度稳定在6000 m左右，起伏不超过20 m，航线经过区域为沙漠地带，地面起伏不大.由EGM2008计算该区域的重力扰动分布如图 4.

采用直接求差法、状态模型法计算的重力扰动结果与EGM2008模型计算的该航线上的重力扰动如图 5，从图中我们可以看出，利用状态模型法计算出的重力扰动和EGM2008模型的趋势是一致的，自适应状态模型法又较状态模型法的趋势更为平滑.直接求差法的结果与整体趋势偏离程度较大，可能是受POS性能的影响(与商用重力仪有一定差距)，计算出的姿态角的精度不是足够高，仍对最后的重力扰动结果带来了影响.表 2为解算结果与EGM2008模型的较差统计表，从表中我们可以看出，利用状态模型法求出的重力扰动与EGM2008模型的差值中误差为8.27 mGal，利用自适应状态模型法的差值中误差为5.43 mGal，接近GB/T17944-2000《加密重力测量规范》困难地区10 mGal的限差要求.

3.2 算例2

算例2为对新疆某区平均航高5000 m、航迹140 km、航速分别为145、135、137 m·s^-1的三条航线进行航空重力解算，飞行轨迹如图 6.飞机在飞行过程中受到的大气扰动较小，三条航线的行高起伏不大.航线经过的区域北部为绿洲，南部为平坦的沙漠，自北向南地面高程逐渐增加，整体无明显起伏.

三条航线的解算结果和EGM2008重力扰动进行作差比较如表 3，从表中可以看出，三条航线差值的中误差都在10 mGal左右，与算例1的结论是一致的.

4 结论与展望

(1)利用状态模型法可以有效地获取高精度、高分辨率的空中重力扰动成果，引入自适应因子以后，结果的精度也得到了一定程度的改善；与EGM2008全球重力场模型对比的差值中误差在10 mGal左右，接近国家在困难地区重力测量的精度要求.

(2)由于采用的惯性系统精度不是很高，未能完全满足直接求差法的要求.从这一点也可以看出状态模型法对设备的性能的要求要低一些.

(3)受缺少地面重力资料的限制，处理结果仅与开放性的重力场模型进行对比验证，在以后的研究中还需要在含有地面重力资料或高程异常的区域进行航空重力测量，来进一步验证本文提出的方法的有效性.

(4)利用航空摄影测量任务的POS数据进行航空重力测量分析，提高了数据的利用率，扩展了航空重力测量的方法，在一定程度上节约了航空重力测量的成本.

(5)基于状态模型法的航空重力测量方法同时还是一种矢量重力测量方法，目前国内还没有成熟的航空矢量重力测量系统，在航空矢量重力数据处理方面也尚处于起步阶段(宁津生，2014)，充分吸收国内外现有研究成果，对我国的重力测量事业有非常重要的意义.