理论模型分析卡尔曼滤波在航空全张量磁力梯度测量中的应用效果

孟庆奎,周坚鑫,舒晴,高维,徐光晶,王晨阳

(1.自然资源部 航空地球物理与遥感地质重点实验室，北京 100083; 2.中国自然资源航空物探遥感中心，北京 100083)

0 引言

实验室条件下，基于超导量子干涉器(superconducting quantum interference device，SQUID)的全张量磁力梯度测量系统的灵敏度可以达到几个fT的量级，能实现对微弱磁性目标体的定位与识别[1-2]，为目标探测提供有效手段。21世纪以来，航空全张量磁力梯度测量系统研发及配套数据处理解释理论是国内外相关领域学者研究的热点，德国、美国、澳大利亚等发达国家对全张量磁力梯度测量技术开展了大量探索研发工作。目前，国际上仅德国光学分子研究中心拥有实用化水平的航空全张量磁力梯度测量系统[3]，且该技术属于西方国家限制出口的尖端技术，公开的文献资料也非常有限。面对技术封锁，我国自力更生，在“国家重大科研装备研制项目”和“国家863计划”支持下，中国科学院上海微系统与信息技术研究所及吉林大学研发团队分别研制出基于低温和高温SQUID的航空全张量磁力梯度测量系统，研发团队攻克了串扰、平台震动监测、磁张量补偿等多项技术难关[4-9]，取得了长足进步，但离实用化还存在一定的差距。我国学者在全张量磁力梯度数据处理[10-11]、正反演[12-15]、综合应用[16-19]等方面也取得了系列成果。

航空全张量磁力梯度测量数据中包含非常复杂的运动噪声，在频谱图中从低频至高频均有分布，并以白噪声为主，如何有效压制运动噪声是一个较大的挑战[3]。Stolz R总结了运动噪声来源，包括仪器本身几何位置变化、外部环境噪声、飞机发动机噪声以及仪器本证噪声[3]，航磁软补偿技术可以去除运动载体在地磁场中产生的恒定磁场、感应磁场和涡流磁场等干扰场[20]，对于白噪声却是无能为力。传统的数字滤波只能滤除指定频段的噪声，对于混叠在全张量磁力梯度有用信号中的噪声不能有效分离。卡尔曼滤波技术在航空重力测量中得到了应用，实践表明，该方法可以有效提取淹没在载体高频振动噪声中的有用信号[21-23]，将其推广应用到航空全张量磁力梯度数据处理中十分必要。张兴东将卡尔曼滤波应用到全张量磁力梯度数据处理中[9]，但是没有明确状态方程和观测方程的搭建方式，而不合理的状态方程会直接影响到卡尔曼滤波效果。笔者借鉴时间序列分析思想，基于最小差分方差建立状态方程，并将卡尔曼滤波与数字滤波相结合，通过模型计算验证了方法的有效性，全区噪声衰减因子优于0.92，即能够去除92%以上噪声成分，全区均方误差优于10 pT/m，满足航空全张量磁力梯度测量精度要求。卡尔曼滤波是一种快速、高效的递推算法，利用实时测量信息对估计值进行最优估计，可直接应用于航空全张量磁力梯度测量实时数据处理。

1 卡尔曼滤波方程建立

1.1 卡尔曼滤波基本理论

滤波估计经历了最小二乘、温纳滤波和卡尔曼滤波的发展而不断完善。1960年，卡尔曼提出离散系统卡尔曼滤波，次年与布西合作，将这种滤波方法推广到连续时间系统中，从而形成了卡尔曼滤波设计理论[24]。卡尔曼滤波是一种时间域滤波方法，采用状态空间方法描述系统，算法采用递推形式，数据存储量小，不仅可以处理平稳随机过程，也可以处理多维和非平稳随机过程。

卡尔曼滤波采用如下空间模型描述动态系统：

X(k)=ΦX(k-1)+ΓW(k-1)，

(1)

Y(k)=HX(k)+V(k)，

(2)

式中：k为离散时间，系统在时刻k的状态为X(k)；Y(k)为对应状态的观测信号；W(k)为输入的白噪声；V(k)为观测噪声；Φ为状态转移矩阵；Γ为噪声驱动矩阵；H为观测矩阵。式(1)为状态方程，式(2)为观测方程。

(3)

P(k,k-1)=ΦP(k-1)ΦT+ΓQΓT，

(4)

K(k)=P(k,k-1)HT[HP(k,k-1)HT+R]-1，

(5)

(6)

P(k)=[I-K(k)H]P(k,k-1)，

(7)

1.2 全张量磁力梯度数据卡尔曼滤波方程建立

全张量磁力梯度仪在每个采样点测量6个梯度值，定义卡尔曼滤波器的状态序列为：

X=(Bxx,Bxy,Bxz,Byy,Byz,Bzz)T，

(8)

各个测量时刻的测量序列为：

Z=(Bxx,Bxy,Bxz,Byy,Byz,Bzz)T，

(9)

因为测量序列与状态序列变量相同，所以可以容易的定义观测矩阵为：

H=diag(1,1,1,1,1,1)，

(10)

式中：diag表示对角矩阵。状态转移矩阵Φ的定义是一个难点，因全张量磁力梯度仪测量的数据是随时间变化的变量，故可借鉴时间序列分析思想，计算包含白噪声的梯度测量值的各阶次差分方差，方差最小的阶次差分的残差可视为白噪声，根据该阶次差分方程可构建状态方程。经大量模型计算，一阶差分方差最小，故可定义状态转移矩阵为：

Φ=diag(1,1,1,1,1,1)。

(11)

本文中卡尔曼滤波进行迭代的本质是利用两个正态分布的融合仍是正态分布的特性而进行的，所指的两个正态分布分别为系统输入噪声和测量噪声的正态分布，系统输入噪声是指我们建立的系统状态方程本身存在的系统误差，测量噪声是指测量仪器精度、外界温度变化等因素产生的测量误差。融合后的正态分布所选择的权值是指卡尔曼增益矩阵，其计算公式可参见式(5)。

2 模型算例分析

2.1 球体模型正演

自然界中存在诸多尺度有限的磁性地质体，当其中心埋深远远大于其直径时，地质体所引起的磁异常特征与球体磁场特征基本吻合，所以研究球体的磁场特征具有一定的实际意义。根据重磁位泊松公式，由引力位计算磁场三分量，进而对磁场三分量在3个不同方向(X、Y、Z)求导便可得到球体的全张量梯度。

模型参数：测区范围2 000 m×2 000 m，测量间距5 m，球心位置(1 000 m，1 000 m，550 m)，球体半径500 m，总磁化强度倾角65°，偏角11°，总磁化强度0.5 A/m。正演模拟结果见图1所示。

2.2 卡尔曼滤波效果

在上述正演结果中加入高斯白噪声V(k)，视为观测噪声，V(k)均值为零，标准差为0.032 nT/m。加入噪声的球体模型全张量磁力梯度视为测量值，如图2所示，可见测量结果全区域都受到高斯白噪声影响，且较严重。遵循以上卡尔曼滤波状态方程和测量方程，建立递推公式，对测量结果进行卡尔曼滤波。在具体操作中，为模拟航空全张量磁力梯度实时测量与实时处理，从第1条测线(即剖面Y=0)的第1个测点开始，随着数据的采集，时间的推移，逐点进行卡尔曼滤波处理，直到全部数据采集完毕，同时完成滤波处理，体现了滤波技术的实时性。在不改变异常形态的前提下，为得到更加光滑的滤波效果，卡尔曼滤波后进行小尺度中值滤波，最终的滤波结果见图3～6所示。

为了定量评价本文所提出方法的滤波效果，进一步评价卡尔曼滤波的降噪性能，引入Pajot等[25]定义的噪声衰减因子β来估计噪声衰减量：

β=[var(Bzs-B)-var(Bkf-B)]/var(Bzs-B)，

(12)

式中:Bzs表示含噪声测量值，B表示理论值，Bkf表示卡尔曼滤波后的值。计算6个梯度张量的全区噪声衰减因子，计算结果见图6a所示，可见β值优于0.92，即能够去除92%以上噪声成分。同时，计算6个梯度张量的全区均方误差，结果见图6b所示，均优于10 pT/m。另外，计算了典型剖面(Y=500 m)上6个梯度张量的噪声衰减因子(0.921 2～0.966 8)和均方误差(6.2～8.5 pT/m)，具体结果见表1所示。以上定量评价结果表明，本文所提方法的滤波效果能够满足当前航空全张量磁力梯度测量精度要求。

图1 球体模型全张量磁力梯度Fig.1 Full tensor magnetic gradient of sphere model

图2 加噪球体模型全张量磁力梯度Fig.2 Full tensor magnetic gradient of noisy sphere model

图3 滤波后全张量磁力梯度立体图Fig.3 Filtered full tensor magnetic gradient of sphere model

图4 Bxx、Byy和Bzz滤波前后对比(剖面Y=500 m)Fig.4 Comparison of Bxx ，Byyand Bzz with before and after filtering(profile Y=500 m)

表1 典型剖面(Y=500 m)卡尔曼滤波效果定量统计

图5 Bxy和Bxz和Byz滤波前后对比(剖面Y=500 m)Fig.5 Comparison of Bxy ，Bxz and Byzwith before and after filtering (profile Y=500 m)

图6 全区噪声衰减因子及全区均方误差Fig.6 Full area noise-reduction factor and mean square error

3 结论

1) 本文通过理论模型验证了卡尔曼滤波方法在航空全张量梯度测量数据中的有效性，结果显示全区噪声衰减因子优于0.92，全区均方误差优于10 pT/m。

2) 由于因卡尔曼滤波对数据的使用率低及噪声随机特征影响，单纯应用卡尔曼滤波不能得到平滑的预测结果，需要串联其他适合的平滑滤波方法，达到更好的滤波效果。

3) 卡尔曼滤波可以利用实时测量信息对估计值进行最优估计，具备快速、高效的递推特性，建议应用于航空全张量磁力梯度实时测量数据处理，也可为其他航空地球物理实时测量数据处理提供参考。