1 引言

航空重力测量是以飞机为载体，测量近地空中重力异常的重力测量技术，它综合运用了航空重力仪，全球定位系统GPS(Global Position System)，惯性导航系统INS(Inertial Navigation System)，无线电等技术设备(张昌达，2005；熊盛青，2009a；王静波等，2009，2010).航空重力测量数据经过后处理，便可得到所需的航空重力异常，这个数据后处理的过程称为航空重力异常解算(简称为航重解算).近年来，差分GPS技术的进步使得航空重力测量精度大幅度提高，解算所得重力异常精度大部分可达到1~2 mGal(1 mGal=10^-5 m·s^-2)，最高的已达0.5 mGal(郭志宏等，2008；熊盛青，2009b; 熊盛青等，2010；Cai et al.，2013).

航空重力测量是在高速飞行运动中进行的，测量数据必然会受到包括飞机运动产生的扰动加速度的干扰，干扰噪声的幅度最高可达到10⁶ mGal以上(孙中苗，2004)，远远高于几十到几百毫伽的重力异常信号.航空重力测量数据的噪声主要体现为高频干扰，依据测量数据的频域特点设计低通滤波器可消除噪声影响提取所需的重力异常信号.目前，能够实现线性相位的FIR(Finite Impulse Response)滤波器是航重解算普遍采用的低通滤波器，利用FIR低通滤波器可以获得精度约1 mGal左右的航空重力异常(郭志宏等，2007，2011).FIR滤波器的显著优点是具有成熟的设计和实现方法，但由于重力信号与噪声在频谱上没有明显的分界线，为获得高精度的解算输出，需设计高阶的低通滤波器以分离重力异常信号与噪声(蔡劭琨等，2010；Cai et al.，2010；罗峰等，2012)，这加重了边界效应对解算结果的影响.

部分学者致力于利用卡尔曼滤波器代替数字低通滤波器实现航空重力异常的解算(Bolotin et al.，2003；Bolotin and Popelensky，2007；Bolotin and Golovan，2013).卡尔曼滤波器通常是采用状态空间模型在时间域对系统的运动状态进行描述，因此卡尔曼滤波器是基于系统运动模型的滤波方法.状态空间模型的应用将不可观测的变量(状态变量)并入可观测的模型，并与其一起处理得到相应的估计输出，所以卡尔曼滤波可实现多维信号的并行处理(王静波，2010；王静波等，2011).在进行航重解算时，还能够对航空重力仪校准参数进行实时估计(Zhou and Cai，2013；Zheng et al.，2014).此外，航重解算是在线下进行的，可对卡尔曼滤波的估计结果进行平滑处理，获得高精度的航空重力异常输出(潘炎冰，2012；李瑞，2014).但是，卡尔曼滤波能够获得准确解算输出的重要前提，是能够掌握系统噪声和量测噪声的统计信息，而这一前提往往是难以满足的.解决这一问题的办法就是用自适应卡尔曼滤波器(AKF，Adaptive Kalman filter)代替传统的标准卡尔曼滤波器(SK，Standard Kalman filter)进行航重解算.基于新息协方差估计的自适应卡尔曼滤波器(IAE，Innovation-based Adaptive Estimation)已证明可以提高噪声统计特性不明时的估计精度(Mohamed，1999；Mohamed and Schwarz，1999)，IAE的优点在于实现方法简单且能够对量测噪声协方差进行实时估计，但其新息协方差估计器的滑动采样窗口宽度缺乏明确的确定方法.

航重解算中，卡尔曼滤波所需的量测量由GPS测定的位置或速度数据提供.由于测量过程中飞机运动状态变化剧烈且具有不规律性，量测噪声参数难以准确预设，因此这里着重讨论利用IAE滤波器实现量测噪声自适应的航重解算方法.本文将依据航空重力测量的基本原理，建立航空重力异常解算的卡尔曼滤波模型，设计适用于G型航空重力测量系统的IAE滤波器.另外，针对IAE滤波器的滑动窗口宽度难以准确确定的问题，采用多个不同滑动窗口宽度的新息协方差估计器，通过对多个新息协方差估计的加权平均，获得滑动窗口宽度优化的改进IAE滤波器.最后，将IAE滤波器和改进的IAE滤波器用于实测数据的航重解算试验.

2 新息自适应卡尔曼滤波器及改进2.1 新息自适应卡尔曼滤波器

设在时刻k，离散线性系统的状态方程和量测方程可分别表示为

X_k为k时刻的状态变量；Φ_k，k－1表示由k-1时刻到k时刻的一步转移矩阵；Γ_k－1为系统噪声驱动阵；Z_k为时刻k的量测量；H_k为量测阵；W_k－1为系统激励噪声序列；V_k为量测噪声序列.Q_k是系统噪声的协方差阵；R_k为量测噪声协方差阵，则标准卡尔曼滤波的算法为(秦永元等，2012)

状态一步预测：

一步预测误差：

滤波增益：

状态估计：

估计均方误差：

新息(Innovation)定义为k时刻滤波器的量测估计值与实际量测值之差，它表示由量测量带来的、滤波器无法预测的那部分信息，即

新息序列的理论协方差为C_rk，则有

当系统噪声和量测噪声服从不相关的高斯分布且新息序列满足各态遍历性时，可根据新息序列的滑动平均获得新息协方差的极大似然估计，有(Mohamed and Schwarz，1999；岳晓奎和袁建平，2005；卞鸿巍等，2006)

上式中，j₀=k-N+1，N为滑动采样窗口的宽度.用 代替式(9)中的C_rk得到基于

新息协方差估计的滤波增益计算公式为

则量测噪声协方差阵的自适应估计为

2.2 滑动采样窗口宽度的优化

新息自适应卡尔曼滤波器根据新息协方差的滑动平均估计实现量测噪声的自适应估计.由式(10)可知，滑动窗口宽度N影响着 的估计结果，但滑动窗口宽度的选取并没有明确的确定方法，通常依赖设计者的经验进行选取.在实际应用中，N不宜选取的过小，否则不能得到的准确估计；N也不宜选取的过大，否则会失去新息的动态特性，致使新息协方差的估计不能反映系统运动状态的变化.

实际应用中，如建立的系统模型不准确，不能真实反映物理过程时，易导致模型与获得的量测值不匹配，将会出现滤波发散，这时需通过增大当前时刻量测信息在状态估计中的加权系数来抑制滤波发散(秦永元等，2012).

将以上两种情况进行综合考虑，由式(11)可知，增益与新协方差估计成反比关系，即越小，当前时刻的滤波增益越高.所以较小的新息协方差有助于增大量测信息的加权系数，使得滤波器可以更多地利用量测信息对一步预测进行修正.因此，当建模粗糙或者系统运动状态不稳定时，可为IAE滤波器设计多个不同滑动采样窗口宽度的新息协方差估计器，通过将各新息协方差估计器的估计结果进行加权平均，获得优化的新息协方差估计.在加权平均的过程中，为了保证滤波器能够较多地利用量测信息进行状态估计，为较小的新息协方差赋予较大的权重.

设M个新息协方差估计器的采样窗口分别为N₁，N₂，…，N_M，其对应的新息协方差估计分别为 ，则采样窗口宽度优化的新息协方差估计的逆为

average_ 为加权平均后的新息协方差估计，w_i为对应新息协方差估计的加权系数(覃方君等，2009)，有

tr[·]表示取矩阵的迹的运算符号； 表示第i个新息协方差估计器的估计，则式(11)变为

可见，相比于普通的IAE滤波器，改进的IAE滤波器配备有多个不同滑动采样窗口宽度的新息协方差估计器，通过对不同采样窗口的新息协方差估计进行加权平均，允许新息协方差估计器的滑动采样窗口在一定范围内动态变化，一方面可以减弱不当的滑动采样窗口造成的解算误差，另一方面也可以提高滤波器工作的稳定性，其估计精度劣于取最优滑动采样窗口的情况，但优于其他滑动采样窗口宽度的估计结果.多个新息方差估计器的结构如图 1所示.

2.3 RTS(Rauch-Tung-Striebel)平滑器

RTS平滑是一种固定区间平滑方法，它利用整个时间区间上的量测量对区间内每个时刻的状态进行估计，而滤波只利用当前时刻以前的量测量作最优估计，因此平滑的精度通常高于滤波的精度.RTS平滑的处理过程是对于给定的在t₁到t_n时间区间上的量测量Z₁，Z₂，…，Z_n，先按顺序作正向的卡尔曼滤波处理，获得并顺序记录RTS平滑器所需的 当k= n时，逆序执行RTS平滑算法(秦永元等，2012)：

为平滑后状态估计输出；K_k ^sm为平滑增益；P_k ^sm为平滑误差.

3 航空重力异常解算的卡尔曼滤波模型

G型航空重力测量系统是三轴稳定平台型测量系统，它包括一个具有水平稳定作用的惯性平台(INS)，一个重力传感器(航空重力仪)，以及机载GPS接收器(移动站)和地面GPS接收器(接收站).本部分将根据航空重力仪运动方程及GPS测量方程构建航重解算的卡尔曼滤波模型.

3.1 航空重力测量的基本原理

航空重力测量的基本原理是从航空重力仪测得的比力数据(飞机运动加速度与重力加速度之和)中扣除飞机运动加速度的影响，获得所需要的重力加速度(熊盛青等，2010).根据牛顿第二定律，航空重力仪(垂直加速度计)传感器探测质量(SM，Sensitive Mass)在地理垂直方向上的动力学方程可表示为(Bolotin，2009)

为重力仪在地理垂直方向上的加速度(高程h的二次导数)；f₃为作用在SM上的比力的垂直分量；ΔG_E为厄特缶斯改正项；G₀(φ，h)为正常场和高度改正项，φ为飞机所处地理纬度；Δg是沿飞行测线的自由空气重力异常(即航空重力异常).

厄特缶斯改正是为了补偿因地球自转产生的离心力和飞机运动产生的离心力附加在重力仪测量质点造成的测量误差，可计算为

u_E为地球自转速度；φ为飞机所处地理纬度；v_E 、v_N 分别是飞机的东向速度和北向速度；R_E 、R_N 分别是地理东向和北向的曲率半径.

重力正常场和高度改正项G₀(φ，h)可由Helmert公式近似计算(Bolotin，2009)：

式中，G=9.7803 m·s^-2，β=0.005302，

β₁=－0.000007，β₀=－0.00014 m·s^-2，

ω₀=0.00123 s^-1.

不难看出，为获得式(19)中的重力异常Δg，需要对航空重力仪测量数据进行一系列的改正.其中，一部分改正项有严密的解析公式，可以直接计算得到，如厄特缶斯改正、正常场和高度改正，它们统称为解析类改正；而另一类改正无法直接用解析公式来计算，需借助滤波方法来实现，即垂直加速度改正.因此，利用卡尔曼滤波实现重力异常解算实际上是在消除噪声干扰的同时实现垂直加速度改正的过程.

3.2 卡尔曼滤波模型的建立

将航空重力仪SM的实时高程，作为卡尔曼滤波的量测方程，写为

h′为重力仪SM所处的实时高度；h为SM的实际高度；δh为测量误差.

航空重力仪及两个水平加速度计的测量方程可表示为(Bolotin，2009)

式(23)中，f′₁、f′₂分别表示水平加速度计测得的东向和北向的运动加速度；f_z1、f_z2为东向和北向实际运动加速度；δf₁、δf₂分别为东向和北向运动加速度的测量误差；f′₃为航空重力仪的测量值，δf₃为测量误差；f_z3为作用在SM的比力在重力仪敏感轴方向上的投影分量；α₁，α₂是INS平台偏离水平状态的误差角；KF₁，KF₂是重力仪安装在INS平台上的角度误差；KF₃是重力仪的尺度因子；TF₃是系统内部的时间延迟.KF₁，KF₂，KF₃，TF₃通称为重力仪校准参数，它们在一次测量过程中通常是固定不变的(Bolotin and Golovan，2013)，可利用卡尔曼滤波对它们进行估计.

将式(23)代入式(19)联立并整理得

其中

f_Σ为已知的作用在SM上的合力产生的加速度总和；q_Σ为系统噪声项.为实现航重解算，将重力异常看作沿测线分布的随机过程，并用滤波器表示，有

将式(24)和式(26)合并作为卡尔曼滤波的系统方程，式(22)为卡尔曼滤波的量测方程，并用状态空间形式表示为

其中，X_f为状态向量，重力异常Δg为其分量，即

A_f为动力学矩阵；B_f为系统噪声驱动矩阵；C_f为控制项驱动矩阵；q_f和δh分别为系统噪声和量测噪声项；u_f为系统的固定控制项.将式(27)和(28)进行离散化，再利用设计的IAE滤波器和RTS平滑器进行处理，即可实现重力异常的估计.

4 航空重力异常解算试验4.1 测量基本情况

我们在2007年利用G型航空重力测量系统在某地进行了实地测量.实验所用数据为第五架次飞行测量数据，飞行时间为4 h.此次测量是为了验证系统工作的稳定性及重复测量结果的一致性.图 2为测线位置示意图，图中横轴为相对经度，纵轴为相对纬度.由图可知，此次测量共有6条东西方向的飞行测线，其中1号、2号、3号及6号测线为四条短测线，它们是沿同一位置往返飞行4次的重复测线；4号测线和5号测线为两条长飞行测线，其中5号测线相较其他测线有约0.002°的北向偏移.测量选取在气象条件较好的夜间进行，以减小气流对测量结果的影响.测量中，设计飞行高度为400 m，实际飞行测量高度在设计飞行高度的±15 m范围内波动，图 3为飞行测量高度变化示意图.此外，飞行测量中的设计飞行速度为60 m·s^-1.

剔除实测数据中的非测线数据，并进行解析类改正后，所得结果(未滤波)如图 4所示.不难看出，航空重力测量信号受噪声干扰非常严重，几十到几百毫伽的重力异常完全被数万毫伽的噪声所淹没.

关于解算试验有以下两点需要说明，第一，由于暂无法获得实际地面重力测量数据，试验过程中将航空重力测量系统配套软件解算结果作为解算参考标准，将本试验的解算结果与参考标准进行对比研究；第二，由于RTS平滑后的重力异常精度高于只使用卡尔曼滤波器估计的重力异常精度(潘炎冰，2012；李瑞，2014)，因此后文解算结果如无特别说明均是指平滑处理后的重力异常.

4.2 不同预设量测噪声的航空重力异常解算试验

将SK解算的重力异常与参考标准进行对比，通过反复调整量测噪声和系统噪声，获得与解算参考标准误差最小时的量测噪声和系统噪声分别为R₀=(0.0153)²和Q₀=diag[0，(3.5×10^-5)²，0，0，0，0，0，(0.85×10^-7)²]，diag[·]为取对角线元素符号.实验中SK和IAE的系统噪声设定为Q₀，初始量测噪声R分别预设为R/R₀=0.01、0.25、25、100进行解算.同时，为确保新息协方差估计的准确性，IAE的滑动采样窗口宽度暂设为80.

SK和IAE在不同预设量测噪声时，解算的重力异常如图 5所示.为方便分析，将解算所得6条飞行测线的重力异常按时间顺序连接组成一条曲线，红色曲线为解算参考标准，蓝色曲线为SK解算的重力异常曲线，黑色曲线为IAE解算的重力异常曲线.由图可见，SK解算的重力异常曲线随着预设R的变化出现较大的波动，而且预设的R与R₀的偏差越大，SK解算重力异常的误差也越高；而在同样的量测噪声条件下，IAE解算的重力异常曲线基本保持稳定，即IAE可降低因预设量测噪声不当引起的解算误差.

统计SK和IAE解算输出与参考标准的误差，如表 1所示.可以看到，在预设量测噪声R=0.01R₀时，SK解算结果的均方误差高达3.423367 mGal；当预设量测噪声R=0.25R₀时，SK解算重力异常与参考标准的均方误差为1.198765 mGal.在IAE解算误差一栏中，不同预设量测噪声条件下，IAE解算重力异常的最大均方误差为0.923542 mGal，最小均方误差为0.887827 mGal，两者的差距仅为0.035715 mGal，而且IAE解算结果与解算参考标准的误差都保持在1 mGal以下.所以，IAE滤波器在解算过程中对量测噪声协方差进行实时的估计和修正，减弱了解算精度对先验量测噪声统计信息的依赖.因此，相比基于SK的航重解算方法，利用IAE进行航重解算更具有普遍性的意义.

4.3 不同滑动窗口的航空重力异常解算试验

由IAE算法步骤可知，滑动采样窗口宽度N直接关系到新息协方差估计的结果，不当的滑动采样窗口也会导致IAE解算精度的下降.本部分将利用不同滑动窗口的IAE滤波器和改进IAE滤波器进行解算试验，验证改进IAE滤波器的优化效果，并讨论N对IAE解算重力异常精度的影响.试验中，IAE滤波器的滑动采样窗口取N=40、80、120、160、200分别进行解算；改进IAE滤波器的新息协方差估计器设定为5个，各新息协方差估计器的采样窗口分别为N₁=40，N₂=80，N₃=120，N₄=160，N₅=200.

图 6为不同滑动窗口的IAE滤波器估计的量测噪声协方差.由图可见，在测线上进行飞行测量时，飞机的运动状态平稳，此时量测噪声协方差估计曲线变化平缓；当飞机由一条测线转飞向另一条测线时(A、B、C、D、E指示的部分)，飞机的运动状态变化剧烈，量测噪声协方差估计曲线出现较大的跳动，而且曲线的跳动频率和幅度随着滑动窗口宽度的增加而减小，即滑动窗口越宽，IAE对量测噪声的估计越平稳，反之，IAE估计的量测噪声曲线跳变越剧烈.

统计不同滑动窗口的IAE与改进IAE的解算误差，如表 2所示.不难看出，在N=80时，IAE的解算输出与解算参考标准的误差最小，此时的滑动采样窗口宽度兼顾了系统的运动状态变化和新息协方差估计的准确性；而在其他采样窗口宽度条件下，IAE的解算结果误差均高于N=80的情况.这也印证了前述的滑动窗口宽度选取的不宜太大也不宜太小.

注意到利用改进IAE解算的重力异常与参考标准的均方误差在1.021047 mGal，解算精度劣于N=80时IAE的解算结果，但优于其他的滑动窗口宽度的情况.因此，改进的IAE可以减小由于滑动采样窗口宽度选取不当引起的解算误差.

5 结语

本文利用新息自适应卡尔曼滤波器及改进的新息自适应卡尔曼滤波器实现航空重力异常解算，并得出如下结论：

(1)新息自适应卡尔曼滤波器的应用，可以降低量测噪声协方差选取不准确造成的解算误差，减弱了解算精度对量测噪声先验统计信息的依赖性；

(2)对新息自适应卡尔曼滤波器进行改进，即在新息自适应卡尔曼滤波器中设计多个不同采样窗口的新息协方差估计器，实现新息协方差估计器的滑动采样窗口优化.解算结果表明，改进的新息自适应卡尔曼滤波的应用能够减小由滑动采样窗口选取不当造成的解算误差.

为获得精度更高、更有普遍意义的航空重力异常解算方法，建议进一步加强对自适应卡尔曼滤波器设计方法的研究.同时，还应继续加强对航空重力测量系统的研究，进一步完善航空重力异常解算的卡尔曼滤波模型，提高解算精度.

致谢 感谢中国国土资源航空物探遥感中心的熊盛青、周坚鑫、郭志宏和周锡华对本研究的大力支持和帮助.