CN104062687A - 一种空地一体的地磁场联合观测方法及系统 - Google Patents

Abstract

本发明公开了一种空地一体的地磁场联合观测方法及系统，所述方法包括：在地磁场测量任务开展前，利用移动地面站采集的地磁场信息，确定局部地磁场模型的模型参数；在执行地磁场测量任务时，利用机载部分在飞行航路上采集的地磁场信息和所述模型参数，得到沿飞行航路上的动态模型；利用机载部分获取的载机按预定姿态飞行时的姿态测量误差，得到沿飞行航路上的测量误差模型；利用所述动态模型和所述测量误差模型，得到飞行航路上采集的地磁场信息在地理坐标系中的估计。本发明能够提高在地理坐标系中表示的地磁场向量的测量精度，并能够进行地磁场异常检测。

Description

一种空地一体的地磁场联合观测方法及系统

技术领域

本发明涉及地磁场测量技术，特别涉及一种空地一体的地磁场联合观测方法及相关的系统。

背景技术

地震、海啸、泥石流和火山喷发对人类造成巨大威胁，据报导孕震期间，本地地磁场合平静太阳日日变化会偏离正常值，这要求对地磁场连续的监测，通过分布式网络，由不同平台组成，磁场向量的获得，在机载平台上，需要在统一的地理坐标系中表示，然而，姿态测量包含噪声，考虑到外围地磁场数值范围大，从两极的60000nT到赤道的20000nT。从机体坐标系到地理坐标系的变换会引入不可忽略的误差，这会导致测量异常失败。

几十年来，航空地磁测量已经广泛的应用于各种领域如探矿、油气和地热勘探，和反潜作战中的潜艇探测，所有这些都是以平台为中心的探测，然而只有磁异探测要求是网络节点探测。大部分航空磁测只要求磁场总场测量，然而，对于地震监测用于的地磁场监测是不够的。

下面相关资料致力于航空地磁观测和地磁数据处理。

1、Moll等获得的No.535,5313美国专利，使用神经网络计算basement rock的深度。

2、Moll等获得的No.588,4256美国专利，阐述了一个数据处理系统，它使用神经网络来计算basement rock的深度。以及神经网络的设计方法和工作过程。以上两个专利都没有对磁场测量和磁场转换进行论述。

3、Shiells等获得的No.6,021,577美国专利,阐述了如何测量矿井的方位，它通过一个基站和移动站来测量，基站和移动站都由两到三个相互垂直的磁通门磁强计构成。

4、www.Picoenvirotec.com网站中题名为IMPAC Integrated Multi-Parameter AirborneConsole的文章，阐述了一个机载地磁观测系统，它装备有4个铯光泵磁强计和一个三轴磁通门磁强计，用于提供载机姿态和载机磁环境补偿模型，然而它没有涉及三轴磁通门的输出量测值如何转换到地理坐标系中。

载机电磁环境和太阳辐射引起的磁场日变化，是两个影响载机进行精确地磁测量的重要因素，需要对其进行补偿，本发明中假设载机磁环境和太阳日变化已被补偿，因此只需关注于减少不精确姿态测量误差引起的磁场测量误差，而安装误差和设备误差等系统误差假设已被校正。

为了减少坐标转换带来的磁场误差，一种方法是采用三个以上天线非共线安装的双频GPS接收机来提高姿态测量精度，此种方法保证关联在天线上的测量基线长度满足测量精度要求。

引用的参考文献如下：

1.Jankowski,J.,and C.Sucksdoroff,,Guide for Magnetic Measurements andObservatory Practice,International Association of Geomagnetism and Aeronomy,1996,Warsaw.

2.Alldredge L.R.,Rectangular harmonic analysis applied to the geomagneticfield,JRG,1981,86(B4),pp3021～3026.

3.Manoj,C.,and S.Maus,Observation of Magnetic fields Generated by Tsunamis,Eos.,Vol.92,No.2,11January 2011.

4.Tsi,J.Digital Techniques for Wideband Receivers(2^nd Ed),Norwood,MA02062,U.S.A.,Artech House,Inc.,2001.

5.M.Idan,Nonlinear Smoothing Identification Algorithm with Application to DataConsistency Checks,J.Guidance,Control and Dynamics,16(2),April1993,pp.337～345.

6.Working Group,Concise Mathematic Handbook,Shanghai Educational Press,1977(in Chinese).

7.Kalman,R.,Discovery and Invention:The Newtonian Revolution in SystemsTechnology,J.of Guidance,Control and Dynamics,Vol.26,No.6,Nov.2003,pp.833～837.

发明内容

本发明的目的在于提供一种空地一体的地磁场联合观测方法及系统，能更好地解决地磁场观测数据精度不高的问题。

根据本发明的一个方面，提供了一种空地一体的地磁场联合观测方法，包括：

在地磁场测量任务开展前，利用移动地面站采集的地磁场信息，确定局部地磁场模型的模型参数；

在执行地磁场测量任务时，利用机载部分在飞行航路上采集的地磁场信息和所述模型参数，得到沿飞行航路上的动态模型；

利用机载部分获取的载机按预定姿态飞行时的姿态测量误差，得到沿飞行航路上的测量误差模型；

利用所述动态模型和所述测量误差模型，得到飞行航路上采集的地磁场信息在地理坐标系中的估计。

优选地，在待测区域内选取若干个观测点，利用移动地面站采集所述若干个观测点的地磁场信息，并利用所述若干个观测点的地磁场信息，估计所述待测区域内的局部地磁场模型的模型参数，其中，所述地磁场信息包括地磁场总量和向量，所述模型参数包括所述局部地磁场模型的直流分量和多个谐波分量的系数。

优选地，所述若干个观测点均匀分布在一条基线上，所述基线是一条原点在参考点x₀，且与磁北方向重合的长度为L的直线。

优选地，利用所述局部地磁场模型的模型参数，对所述机载部分在飞行航路上采集的地磁场信息进行分解，得到各个方向的直流分量和谐波分量，并将各个方向上的谐波分量进行离散化处理，得到所述飞行航路上的动态模型。

优选地，机载部分通过其定姿定位系统获取姿态测量误差。

优选地，利用所述动态模型和所述测量误差模型，得到在地理坐标系中的卡尔曼滤波器，并通过所述卡尔曼滤波器，得到飞行航路上采集的地磁场信息在地理坐标系中的估计。

优选地，还包括：利用所述卡尔曼滤波器内的数据流，检测地磁场是否发生异常。

根据本发明的另一方面，提供了一种空地一体的地磁场联合观测系统，包括：

移动地面站，用于在地磁场测量任务开展前，利用移动地面站采集的地磁场信息，确定局部地磁场模型的模型参数；

机载部分，用于在执行地磁场测量任务时，利用在飞行航路上采集的地磁场信息和所述模型参数，得到沿飞行航路上的动态模型，利用获取的载机按预定姿态飞行时的姿态测量误差，得到沿飞行航路上的测量误差模型，并利用所述动态模型和所述测量误差模型，得到飞行航路上采集的地磁场信息在地理坐标系中的估计。

优选地，所述机载部分利用所述动态模型和所述测量误差模型，得到在地理坐标系中的卡尔曼滤波器，并通过所述卡尔曼滤波器，得到飞行航路上采集的地磁场信息在地理坐标系中的估计。

优选地，还包括：

磁异常检测器，用于利用所述卡尔曼滤波器内的数据流，检测地磁场是否发生异常。

与现有技术相比较，本发明的有益效果在于：

本发明结合地球物理学、航电集成技术、磁传感器技术、数字信号处理技术、数字信号处理技术等，能够直接面向空天地一体化观测体系中的航空平台地磁场测量，这些量测数据精度高，可供地震预报和矿藏勘探分析。

附图说明

图1是本发明实施例提供的空地一体的地磁场联合观测方法原理框图；

图2是本发明实施例提供的机载部分组成框图；

图3是本发明实施例提供的移动地面站组成框图。

具体实施方式

以下结合附图对本发明的优选实施例进行详细说明，应当理解，以下所说明的优选实施例仅用于说明和解释本发明，并不用于限定本发明。

本发明提供了一种空地一体的地磁场联合观测系统及相应的数据处理方法。

(1)所提出的空地一体的地磁场联合观测系统分为移动地面站(即移动地面基站)和机载部分(即载机观测控制台或载机观测台)两部分。其中，所述移动地面站在地磁场测量任务开展前，利用移动地面站采集的地磁场信息，确定局部地磁场模型的模型参数；所述机载部分用于在执行地磁场测量任务时，利用在飞行航路上采集的地磁场信息和所述模型参数，得到沿飞行航路上的动态模型，利用获取的载机按预定姿态飞行时的姿态测量误差，得到沿飞行航路上的测量误差模型，并利用所述动态模型和所述测量误差模型，得到飞行航路上采集的地磁场信息在地理坐标系中的估计，具体地说，所述机载部分利用所述动态模型和所述测量误差模型，得到在地理坐标系中的卡尔曼滤波器，并通过所述卡尔曼滤波器，得到飞行航路上采集的地磁场信息在地理坐标系中的估计。进一步地，所述系统还包括磁异常检测器，其利用所述卡尔曼滤波器内的数据流，检测地磁场是否发生异常，所述磁异常检测器可以设置在移动地面站中，也可以设置在机载部分。

所述载机观测台包括一台铯光泵磁强计、一台三轴磁通门磁强计、一套机载定姿定位子系统(POS)和一个综合任务处理单元。铯光泵磁强计用于测量地磁场总场值，三轴磁通门磁强计用于在机体坐标系统中测量地磁场三分量，它们均安装在无磁尾杆中。机载观测台用于在飞行航线上采集地磁场向量和总量测量数据，采用(4)和(6)给出的模型和算法给出地磁场向量在地理坐标系中的估计，同时利用该算法获得的新息，检测可能出现的地磁异常。

所述地面基站，是可移动的，用于模型建模的数据处理单元和飞行后数据处理，一般安置置放于任务区域的中心。在单轴无磁转台上安装有跟机载观测控制台同样的磁传感器，用于飞行前数据采样，建立本地地磁场模型，采用(2)和(6)给出模型和算法，识别出模型的参数。

(2)用具有直流分量的拟ARMA过程来描述文献[2]提出的局部地磁场模型，并能够辨识局部地磁场模型参数，即在地磁场测量任务展开前，利用地面站测量到的磁场数据进行模型参数辨识。

(3)飞行器根据事先设定的直线航路飞行并采集航路点上的地磁场向量和总量。机载观测控制台上的综合任务处理单元4对这些数据进行处理和存储。

(4)利用(2)得到的局部地磁场模型的参数，将机载传感器得到测量分解为直流分量和谐波分量，将各个方向上的谐波分量用状态方程表示，离散化后，提出飞行路线上一个递推形式的局部地磁场动态模型，它用于估计器设计。假设载机在水平面以固定速度沿着给定的航向运动，考虑姿态传感器测量误差后，给出了在上述传感器配置下测量误差模型(即观测模型)。

(5)为了平滑姿态传感器的噪声，考虑到姿态输出的数据率远高于磁传感器的数据率，还可以给出了一个简单的α-姿态滤波器(式(17))，用于改善姿态测量精度。

(6)设计了一个扩展卡尔曼滤波器(EKF)，利用(2)识别出的地磁场模型参数，(4)给出的动态地磁模型和观测模型，(5)给出的姿态滤波器的输出，对(1)中光泵磁强计输出的测量点处地磁场总量和磁通门磁强计输出的测量点处地磁场在机体坐标系中的三个分量，给出在一个统一的地理坐标系下表达的地磁场强度向量的估计。

(7)确定(6)中滤波器初始条件、系统噪声协方差矩阵和测量噪声协方差矩阵。

图1是本发明实施例提供的空地一体的地磁场联合观测方法原理框图，如图1所示，包括：

步骤S101、在地磁场测量任务开展前，利用移动地面站采集的地磁场信息，确定局部地磁场模型的模型参数。

具体地说，在待测区域内选取若干个观测点，利用移动地面站采集所述若干个观测点的地磁场信息，并利用所述若干个观测点的地磁场信息，估计所述待测区域内的局部地磁场模型的模型参数，其中，所述地磁场信息包括地磁场总量和向量，所述模型参数包括所述局部地磁场模型的直流分量(式(2)中的c)和多个谐波分量的系数(式(3)中的a₁...a_K)。其中，所述若干个观测点均匀分布在一条基线上，所述基线是一条原点在参考点x₀，且与磁北方向重合的长度为L的直线。

步骤S102、在执行地磁场测量任务时，利用机载部分在飞行航路上采集的地磁场信息和所述模型参数，得到沿飞行航路上的动态模型。

利用所述局部地磁场模型的模型参数，对所述机载部分在飞行航路上采集的地磁场信息进行分解，得到各个方向的直流分量和谐波分量，并将各个方向上的谐波分量进行离散化处理，得到所述飞行航路上的动态模型(式(11))。

步骤S103、利用机载部分获取的载机按预定姿态飞行时的姿态测量误差，得到沿飞行航路上的测量误差模型(式(16))。其中，机载部分通过其定姿定位系统获取姿态测量误差。

步骤S104、利用所述动态模型和所述测量误差模型，得到飞行航路上采集的地磁场信息在地理坐标系中的估计(式(18))。

具体地说，利用所述动态模型和所述测量误差模型，得到在地理坐标系中的卡尔曼滤波器，并通过所述卡尔曼滤波器，得到飞行航路上采集的地磁场信息在地理坐标系中的估计。

进一步地，还可以利用所述卡尔曼滤波器内的数据流检测地磁场是否发生异常。

图2是本发明实施例提供的机载部分组成框图，如图2所示，机载观测台任务是，对一个以地面站为中心，约3000km×1000km(纬度×经度方向)的矩形区域进行地磁场测量，由矩形区域地形决定在300m到1000m之间的安全高度区域。

图2示意了机载地磁观测台的构成，三轴磁通门磁强计1和铯光泵磁强计2安装在无磁尾杆的末端。为了避免对飞机原有空气动力学性能造成太大影响，它尽可能的长，同时作为桁梁(它的形变程度取决于它的材料)要求在远端的形变程度小于载机的航向和横滚测量误差。为了减少供电电流感应的磁场的干扰，通过布放在尾杆内的双绞线对三轴磁通门磁强计1和铯光泵磁强计2使用机上直流电源供电，并选择在机舱合适位置的地板上进行接地，接地电阻小于10Ω。

机载定姿定位子系统3，是一个商业货架产品，安装在飞机的设备舱内，用来提供机体坐标系中表示的位置、姿态；采样时间以UTC(universal coordinated time)为标准。

一个机载综合任务处理单元4和监视器，安装在飞机操作舱内，用来处理、存储和显示采集到的磁场数据和载机的飞行参数。机载综合任务处理单元4主要运行它设计有合理的人机界面，可实时操作和查看采集的观测数据。机载综合任务处理单元4担负的主要处理任务包括：1)实现本发明提出的地磁场估计(式(11)～式(19))和磁异常检测(式(20))；2)载机磁环境补偿算法；3)地磁日变补偿算法。

图3是本发明实施例提供的移动地面站组成框图，如图3所示，基站和机载观测台的目的在于地磁变化的观测和数据处理。基站可以独立工作对本地地磁场进行测量，三轴磁通门磁强计5和铯光泵磁强计6被安装在一个无磁转台上，可以手动校准，而且平台具有和POS一致的测角精度，同时一个差分GPS接收机(DGPS)7被安装在平台上，用于测量基站的经纬度和数据采样时间。一个地面数据处理单元8，由工控机和显示器组成，用于离线数据存储、管理和处理功能，例如地磁昼变化监控和数据滤波，参数辨识和局部地磁场模型验证。

一、划定一个在郊区的调查区域，通过下述实验方法，可建立一个半经验的局部地磁场模型，为了避免太阳辐射引起的地磁日变，建议在深夜展开实验。

1、选择一条原点在参考点x₀，长度为L和磁北方向重合的直线作为基线l，依据当地地形从2km到数km。沿着基线用DGPS接收机选择间距为Δl的偶数个观测点，x₁,…,x_M,M≥2。在每个观测点x_i地磁场被测量N次(比如N＝100)，其均值作为处理模型的输入。

${\stackrel{\‾}{B}}_x\left(x_i\right)=\frac{1}{N}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{B}_{\mathrm{xn}}\left(x_{\mathrm{im}}\right)$

${\stackrel{\‾}{B}}_z\left(x_m\right)=\frac{1}{N}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{B}_{\mathrm{zn}}\left(x_m\right)$

式中，B_xn(x_m),B_zn(x_i)分别为三轴磁通门磁强计在一个水平轴垂直轴向的在x_i点的第n个测量值，n＝1,…N,m＝1,…M。

2、地磁场可以被分解为一个常量和一系列谐波分量之和([2])，表达式有同样的形式。

${\stackrel{\‾}{B}}_j\left(x_i\right)=c_j+s_j\left(i\right)+v_j,j=x,y,z$

c_j是待定的常量，v_j是测量噪声，其均值非零，协方差σ_vj由磁传感器噪声和在当前点上的采样次数确定：

${\mathrm{\σ}}_{\mathrm{vj}}=\frac{{\mathrm{\σ}}_{\mathrm{sensorj}}}{\sqrt{N}}$

在每个轴向上，谐波分量被近似为一个AR(auto-regression)模型，系数待定。

s_j(i)＝a_1js_j(i-1)+...+a_Kjs_i(i-K_j)

其中阶数K_j≤M-1是偶数，j＝x,y,z。它的特征多项式是：

$1-a_{1j}{\mathrm{\λ}}^{-1}-...-a_{\mathrm{Kj}}{\mathrm{\λ}}^{-K_j}=0---\left(1\right)$

有特征根n＝1,2,…K_j/2，j＝x,z。在复平面单位圆上，ω_n,j代表着空间的周期属性，Δl是采样区间长度。

为了使表述简洁，用来辨识此模型参数的步骤统一描述如下

u(k)＝c+s(k)+v_n(k) (2)

s(k)＝a₁s(k-1)+...+a_Ks(k-K) (3)

其中，k＝1,…M；M≥K+1；阶数K≥2是一个偶数整数，v_n是一个高斯白噪声，这步骤目的是为了提供一个对常量c，以及(2)、(3)式中的参数估计，此步处理的物理意义在后续的“局部地磁场模型参数辨识”部分阐述。

u＝[u(i),u(i+1),…u(i+M)]^T

p＝[q,a₁,…a_K]^T

是需要辨识的量测向量和参数向量，相对的，上标T代表向量转置，q＝c(1-a₁-…-a_K)，上式写成如下矩阵公式：

u＝A_mp

$A_m=\left[\begin{array}{cccc}1& u(k-1)& ...& u(k-K)\\ 1& u\left(k\right)& ...& u(k+1-K)\\ .& & & .\\ .& ...& ...& .\\ .& & & .\\ 1& u\left(M\right)& ...& u(M-K)\end{array}\right]---\left(4\right)$

可推出如下最小二乘估计：

$\hat{p}={\left(A_m^T{A}_m\right)}^{-1}{A}^{T}u---\left(5\right)$

(2)式中常量c的估计如下表示

$\hat{c}={\hat{p}}_1(1-{\hat{p}}_2-...{\hat{p}}_K)---\left(6\right)$

最后，将j＝x,z,本地磁场测量序列，代入前述的步骤，即返回了需要的对应于特定磁通门测量通道的参数和并利用辨识出的局部地磁模型表示j＝x,y,z。所获得的模型参数被存储在综合任务处理单元4中，用于飞行测量数据的处理。具体地，

$\hat{A}={\hat{c}}_x,\hat{C}={\hat{c}}_z---\left(7\right)$

在地磁观测任务中，载机要求在水平面内以预定航向和速度v_g直线飞行，任务处理单元4通过RS-232串口以10Hz采样处理磁强计和POS系统数据，并且存储到综合任务处理单元4中。有两个处理过程被设计但是没有在此发明中描述，用于补偿磁场环境和地磁日变化。

构造了一个用于在线数据处理的局部地磁动态模型被建立，描述如下。

假设载机以速度v_g在水平面内直线飞行，相对本地磁场磁北方向的航向为对磁通门每个通道的输出，采用后续的“局部地磁场模型参数辨识”部分的记号(A7)，记直流常数三个磁通门通道分别为和c_z.A和c_z由式(7)识别得到。则三个通道输出的谐波分量分别为为和s_z＝B_bz--c_z。

特别是，在原点在载机上，指向和以地面站为原点为地理坐标系相同，方向与载机航向一致时，在磁通门每个通道上，离散化后，有

k＝0,1,2, (8)

其中max(K_x,K_y)≤K，Δt是采样时间。不失一般性，设置x₀＝0，ω＝v,Δl＝v_gΔt，对于j＝x,y，和由式(7)确定(见后续的“局部地磁场模型参数辨识”部分)。

然后定义一个状态向量

具有长度N_s＝3+K_x+K_y+K_z,得到如下状态转移方程

x_k＝Φx_k+B_ww_k,k＝0,1,2 (10)

w向量表示一定长AR模型截断误差Φ_Ns×Ns：

$B_w=\left[\begin{array}{ccc}0& 0& 0\\ & ...& \\ 1& 0& 0\\ 0& 0& 0\\ & ...& \\ 0& 1& 0\\ 0& 0& \\ & ...& \\ 0& 0& 1\end{array}\right]$

B_w是N_s×3的实矩阵，在(3+K_x,1),(3+K_x+K_y,2),(N_s,3)处元素非零。均值为0的系统噪声w，其协方差矩阵Q，被假设为半正定矩阵。

载机以给定的姿态θ₀,γ₀飞行，POS系统具有独立的姿态测量误差量测值由y_k＝[B_g+v_g，B+v_B]^T给出，v_g是矩阵协方差的白噪声，v_B是光泵的量测噪声。

$\begin{array}{c}{R}_{B_g}=E\left\{\mathrm{\Δ}{B}_g\mathrm{\Δ}{B}_g^T\right\}\\ =\mathrm{diag}[{\mathrm{\σ}}_{\mathrm{bx}}^2,{\mathrm{\σ}}_{\mathrm{by}}^2,{\mathrm{\σ}}_{\mathrm{bz}}^2]\\ +E\left\{C(B_{g0},{\mathrm{\φ}}_0,{\mathrm{\θ}}_0,{\mathrm{\λ}}_0)\left[\begin{array}{c}\mathrm{\Δ\φ}\\ \mathrm{\Δ\θ}\\ \mathrm{\Δ\γ}\end{array}\right]\right[\mathrm{\Δ\φ},\mathrm{\Δ\θ},\mathrm{\Δ\γ}\left]C^T(B_{g0},{\mathrm{\φ}}_0,{\mathrm{\θ}}_0,{\mathrm{\λ}}_0)\right\}\end{array}---\left(12\right)$

其中：

and是磁通门每个通道上的量测噪声方差，这些know-how在后续的“噪声协方差估计”部分中进行描述。

地磁场在机体坐标系统进行测量：

$B_{\mathrm{bx}}=x_1-x_{3+K_x}{x}_{3+\mathrm{Kx}}$

$B_{\mathrm{by}}=x_2+x_{3+K_x+K_y}$

$B_{\mathrm{bz}}=x_3+x_{3+K_x+K_y+K_z}$

为从机体坐标系到地理坐标系的姿态转换矩阵(见噪声协方差估计部分)，量测值可如下表示：

$y_k=\left[\begin{array}{c}{A}_b^g\left[\begin{array}{c}{x}_1\left(k\right)-x_{3+K_x}\left(k\right)\\ x_2\left(k\right)+x_{3+K_x}\left(k\right)\\ x_3\left(k\right)+x_{3+K_x+K_y+K_z}\left(k\right)\end{array}\right]+v_g\left(k\right)\\ \sqrt{{(x_1\left(k\right)-x_{3+K_x}\left(k\right))}^2+{(x_2\left(k\right)+x_{3+K_x}\left(k\right))}^2+{(x_3\left(k\right)+x_{3+K_x+K_y+K_z}\left(k\right))}^2}+v_B\left(k\right)\end{array}\right]$

$\begin{array}{c}y\left(k\right)=\left[\begin{array}{c}{H}_1{x}_k\\ \sqrt{{(x_1\left(k\right)-x_{3+K_x}\left(k\right))}^2+{(x_2\left(k\right)+x_{3+K_x}\left(k\right))}^2+{\left(x_{3+K_x+K_y+K_z}\left(k\right)\right)}^2}\end{array}\right]+\left[\begin{array}{c}{v}_{\mathrm{gk}}\\ v_{\mathrm{Bk}}\end{array}\right]\\ \≡h_1\left(x_k\right)+v_{\mathrm{mk}}\end{array}---\left(14\right)$

k＝0,1,2,…，

载机在水平面运动时，

v_mk是量测噪声的协方差向量，R_y如下表示：

$R_y=\left[\begin{array}{cc}{R}_{B_g}& 0\\ 0& R_B\end{array}\right]=\left[\begin{array}{cc}{R}_{B_g}& 0\\ 0& E\left(v_B^2\right)\end{array}\right]---\left(16\right)$

同时假设俯仰和横滚角为0，航向不能被假设为一个已知的量，而需要估计出来，因此，假设地磁场模型是与航向完全解耦的运动模型，相较于磁强计自身噪声，此噪声强度较小，在本发明中，它通过如下简单预处理滤波器估计输出：

是测量航向角，从POS系统3中获取。

扩展卡尔曼滤波器([7])(EKF)，嵌入在综合任务处理单元4中执行，用于估计所需的地磁场。在地理坐标系中表示，EKF如下描述：

${\hat{x}}_{k+1|k}=\mathrm{\Φ}{\hat{x}}_{k|k}$

$P_{k+1|k}=\mathrm{\Φ}{P}_{k|k}{\mathrm{\Φ}}^T+B_{w}Q{B}_w^T$

$K_k=P_{k+1|k}\frac{\∂h_1^T\left({\hat{x}}_{k+1|k}\right)}{\∂}{(\frac{\∂h_1\left({\hat{x}}_{k+1|k}\right)}{\∂x}{P}_{k+1|k}\frac{\∂h_1^T\left({\hat{x}}_{k+1|k}\right)}{\∂x}+R_y)}^{-1}$

${\hat{x}}_{k+1|k+1}={\hat{x}}_{k+1|k}+K_k(y_k-h_1\left({\hat{x}}_{k+1|k}\right))$

$P_{k+1|k+1}=(I-K_k\frac{\∂h_1\left({\hat{x}}_{k+1|k}\right)}{\∂x})P_{k+1|k}---\left(18\right)$

给定观察序列{y_i,i＝0,1,2,…,k+1}，表示x_k+1的估计值，是x_k+1的预测值，y_k＝y(k),K_k是滤波器的增益，P_k+1|k和P_k+1|k+1是状态估计和预测的协方差矩阵，k＝0,1,2,…,线性观测函数h₁(x)采用进行一步预测状态。

初始坐标P₀＝P(0)和初始状态方差矩阵采用如下形式，P(0)＝βI，辨识矩阵乘上一个大数β，初始状态向量猜测有一点复杂，前三个分量组成如下表示：

${\hat{x}}_{3+K_x+1}\left(0\right)=y_1\left(1\right)-{\hat{c}}_y$

${\hat{x}}_2\left(0\right)={\hat{c}}_z$

其中，被基站辨识出。剩余其它初始状态向量靠首个从磁强计量测得到K(max(K_x,K_y,K_z))猜测，如下

${\hat{x}}_{3+1}\left(0\right)={\hat{c}}_x-y_1\left(1\right)$

…，

${\hat{x}}_{3+K_x}\left(0\right)={\hat{c}}_x-y_1\left(K_x\right),$

${\hat{x}}_{3+K_x+1}\left(0\right)=y_2\left(1\right)-{\hat{c}}_y,$

…，

${\hat{x}}_{3+K_x+K_y}\left(0\right)=y_2\left(K_y\right)-{\hat{c}}_y,$

${\hat{x}}_{3+K_x+K_y+1}\left(0\right)=y_3\left(1\right)-{\hat{c}}_z,$

…，

${\hat{x}}_{3+K_x+K_y+K_z}\left(0\right)=y_3\left(K_z\right)-{\hat{c}}_z.$

过程噪声方差矩阵Q被从一个3×3的对角矩阵选出，它的元素不小于0.01nT。

本发明还包括磁异常检测器(MAD)，既可以在综合任务处理单元4中，也可以地面数据处理单元8中。该MAD利用从EKF获得的新息统计构建序列，用于检测可能的磁异常。所述的术语新息被用来描述EKF中的数据流如果模型有效，它们都是无偏的。因此，在每个测量通道中，都可以建立检测模型如下：

H0_j:在j方向没有此异常

H1_j:在j方向有磁异常.

在大量采样N_sample数据量的情况下，例如N_sample>100，对于给定的容许漏报概率ε，即地磁异常发生而探测器没有报告的概率，t分布的检测器具有如下形式：

$|{\stackrel{\‾}{y}}_j/{\widehat{\mathrm{\σ}}}_{j\left(k\right)}|t_{\mathrm{\ϵ}},j=\mathrm{1,2,3}.k=N_{\max },N_{\max }+1,...,---\left(20\right)$

时，认定H0_j成立，该分量所对应的地磁场强度无异常，其中

${\stackrel{\‾}{y}}_j\left(k\right)={\mathrm{\Σ}}_{i=1}^{N_{\mathrm{sample}}}{y}_i(k-i-1)/N_{\mathrm{sample}}$

方差的估计如下：

${\widehat{\mathrm{\σ}}}_j\left(k\right)=\sqrt{{\mathrm{\Σ}}_{i=1}^{N_{\mathrm{sample}}}{(y_i(k-i-1)-{\stackrel{\‾}{y}}_i\left(k\right))}^2/(N_{\mathrm{sample}}-1)}$

给定漏报概率ε,门限t_ε可从表中查询到，可从数学手册，例如[6]的t-分布表中查询到。

当测量时滤波数据、传感器数据、采样时间、载机位置和姿态，以预定义的格式存储在综合任务处理单元4磁盘里的一个文件中。载机执行完任务后，数据文件可从综合任务处理单元4下载，被提交于有关机构供地震分析预测。

二、局部地磁场模型参数辨识

众所周知，现有许多模型用于描述地磁场，无论是全球磁场和局部地磁场，例如IGRF([1])、WMM([1])和CHINAMF([4])。就航空地磁观测而言，所关心的区域并没有那么大，而且这些模型对于地震诱发的地磁异常探测来说还不是足够精确，因此，局部地磁场模型来描述磁场看来是合理的，此模型建立于偶极子地磁场假设，由球形谐波函数组成，并且具有足够的精度，可以用于本发明的局部地磁场建模。在本发明中，由外太空辐射诱发的地磁异常被忽略掉，例如太阳风。

以感兴趣的点为原点建立一个地理坐标系，地磁场的势是Laplacian方程的解，可如下描述[2]：

$\begin{array}{c}U-\mathrm{Ax}-\mathrm{By}-\mathrm{Cz}-\underset{q_1=1}{\overset{N_{\max }-1}{\mathrm{\Σ}}}\underset{i=1}{\overset{q_1}{\mathrm{\Σ}}}\{D_{\mathrm{ij}}\cos \left(\mathrm{ivx}\right)\cos \left(\mathrm{jwy}\right)+E_{\mathrm{ij}}\cos \left(\mathrm{ivx}\right)\sin \left(\mathrm{jwy}\right)\\ +F_{\mathrm{ij}}\sin \left(\mathrm{ivx}\right)\sin \left(\mathrm{jwy}\right)+G_{\mathrm{ij}}\sin \left(\mathrm{ivx}\right)\cos \left(\mathrm{jwy}\right)\}\exp \left(u_{1}z\right)\end{array}---\left(A1\right)$

j＝q-i+1,v＝2π/L_x，w＝2π/L_y，u₁＝[(iv)2+(jw)2]^1/2

其中U是局部地磁场的势函数，L_x和L_x是感兴趣的矩形区域长度，A、B和C是常量，D_ij、E_ij、F_ij和G_ij是系数。

在地理坐标系中，地磁观测量如下

B_g＝-μ₀▽U (A2)

μ₀＝4π×10^-7H/m是空气的电介质常数，▽是梯度算子，在一个上述定义的地理坐标系中，

$B_g=\left[\begin{array}{c}A+\underset{q_1=1}{\overset{N_{\max }-1}{\mathrm{\Σ}}}\underset{i=1}{\overset{q_1}{\mathrm{\Σ}}}\mathrm{iv}\{-D_{\mathrm{ij}}\sin \left(\mathrm{ivx}\right)\cos \left(\mathrm{jwy}\right)-E_{\mathrm{ij}}\sin \left(\mathrm{ivx}\right)\sin \left(\mathrm{jwy}\right)\\ +F_{\mathrm{ij}}\cos \left(\mathrm{ivx}\right)\sin \left(\mathrm{jwy}\right)+G_{\mathrm{ij}}\cos \left(\mathrm{ivx}\right)\cos \left(\mathrm{jwy}\right)\}\exp \left(u_{1}z\right)\\ B+\underset{q_1=1}{\overset{N_{\max }-1}{\mathrm{\Σ}}}\underset{i=1}{\overset{q_1}{\mathrm{\Σ}}}\mathrm{jw}\{D_{\mathrm{ij}}\cos \left(\mathrm{ivx}\right)\sin \left(\mathrm{jwy}\right)-E_{\mathrm{ij}}\cos \left(\mathrm{ivx}\right)\cos \left(\mathrm{jwy}\right)\\ -F_{\mathrm{ij}}\sin \left(\mathrm{ivx}\right)\cos \left(\mathrm{jwy}\right)+G_{\mathrm{ij}}\sin \left(\mathrm{ivx}\right)\sin \left(\mathrm{jwy}\right)\}\exp \left(u_{1}z\right)\\ C+\underset{q_1=1}{\overset{N_{\max }-1}{\mathrm{\Σ}}}\underset{i=1}{\overset{q_1}{\mathrm{\Σ}}}{u}_1\{D_{\mathrm{ij}}\cos \left(\mathrm{ivx}\right)\cos \left(\mathrm{jwy}\right)+E_{\mathrm{ij}}\cos \left(\mathrm{ivx}\right)\sin \left(\mathrm{jwy}\right)\\ +F_{\mathrm{ij}}\sin \left(\mathrm{ivx}\right)\sin \left(\mathrm{jwy}\right)+G_{\mathrm{ij}}\sin \left(\mathrm{ivx}\right)\cos \left(\mathrm{jwy}\right)\}\exp \left(u_{1}z\right)\end{array}\right]---\left(A3\right)$

(A3)适用于任何惯性坐标系，特别的当Ox与基线L重合,Oy轴与Ox垂直时，可简化(A3)如下：

$B_g=\left[\begin{array}{c}A+\underset{q=1}{\overset{N_{\max }-1}{\mathrm{\Σ}}}\underset{i=1}{\overset{q}{\mathrm{\Σ}}}\mathrm{iv}\{-D_{\mathrm{ij}}\sin \left(\mathrm{ivx}\right)-E_{\mathrm{ij}}\cos \left(\mathrm{ivx}\right)\}\exp \left(\mathrm{vz}\right)\\ \underset{q=1}{\overset{N_{\max }-1}{\mathrm{\Σ}}}\underset{i=1}{\overset{q}{\mathrm{\Σ}}}\mathrm{jw}\{E_{\mathrm{ij}}\cos \left(\mathrm{ivx}\right)+F_{\mathrm{ij}}\sin \left(\mathrm{ivx}\right)\}\exp \left(u_{1}z\right)\\ C+\underset{q=1}{\overset{N_{\max }-1}{\mathrm{\Σ}}}\underset{i=1}{\overset{q}{\mathrm{\Σ}}}{u}_1\{D_{\mathrm{ij}}\cos \left(\mathrm{ivx}\right)+G_{\mathrm{ij}}\sin \left(\mathrm{ivx}\right)\}\exp \left(\mathrm{vz}\right)\end{array}\right]---\left(A4\right)$

式(A4)表明了地磁场可以局部的分解为一系列谐波分量的组成，其“频率”和系数未知。针对空间“频率”，阐明两个观点是有用的。

首先利用基本的三角关系，(A4)可以被表示为如下更紧凑的形式：

当z≠0,(A4)中参数D_ij和G_ij可以通过自乘e^ivz，i＝0,1,…2来更新，而谐波分量的幅度和相位全部依赖于基站的高度，这种乘法运算不改变(A5)的频率结构。

其次，对任意地理坐标系Ox’y’z’，其方向φ，根据微分几何学，

$\left[\begin{array}{c}\frac{\∂}{\∂x^{\′}}\\ \frac{\∂}{\∂y^{\′}}\\ \frac{\∂}{\∂z^{\′}}\end{array}\right]=\left[\begin{array}{ccc}\cos \mathrm{\φ}& -\sin \mathrm{\φ}& 0\\ \sin \mathrm{\φ}& \cos \mathrm{\φ}& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}\frac{\∂}{\∂x}\\ \frac{\∂}{\∂y}\\ \frac{\∂}{\∂z}\end{array}\right]---\left(A6\right)$

它进一步阐明了(A4)中基频v是独立于特定的坐标系，在新的方向上保持不变，从而，一旦谐波分量的频率被辨识，式(A5)或A(6)同样的可应用于不同高度的地磁场。

(A4)中求和项数为1+2+…N_max-1＝(N_max-1)N_max/2,其数值可能很大，为了确保问题可解，在范围2～(N_max-1)N_max/2内，在实践中，使用偶数个数据K_j，j＝x,y,z，以限制磁通门每个通道上AR模型的长度。因此，在模型精度和数值易处理上折衷是需要的，虽然这种截断会导致很小的模型误差。

依照常规的处理方法([4])，可知一个具有不同的谐波分量组成信号

具有未知但是不变的相位，是下列方程的解。

s(k)＝a₁s(i-1)+...+a_Ks(k-K)k＝1,…M,M≥K+1 (A7)

该方程的特征多项式

1-a₁λ^-1-...-a_Kλ^-K＝0 (A9)

具有在复平面的单位圆上的根n＝1,2,…K_j/2,j＝x,y,z,，代表了每个谐波分量的“频率”。注意到多项式没有实数根

1-a₁-...-a_K≠0 (A10)

注意，式(A5)的每个分量都不具有(A7)的形式。实际上，它们可以被统一表述为具有未知直流分量，受噪声噪声污染的信号模型，

u(i)＝c+s(i)+v_n(i) (A11)

v_n(i)代表了传感器噪声。特别是，谐波分量s(i)不能从量测值中分离出来。

在信号处理的角度，采用通常足够的大时间来覆盖一个周期，如此直流分量c可以容易的被滤去，结果，大部分采用频率估计的方法都假设信号中没有直流分量。然而，基线的长度是远远小于典型观测区域，其长宽L_x和L_x通常在几百到上千公里的范围([1])，式(A11)中的常量不能被滤除，也不能被分离出来。因此挑战在于，使用基线上的小数目采样点来估计未知常量和系数，这是本发明与已有技术最大的区别。

为了处理这种麻烦的事情，设计了一种处理过程如下。注意到

u(k)＝c+s(k)+v_n(k)

a₁u(k-1)＝a₁c+a₁s(k-1)+a₁v_n(k-1)(A12)

·

·

·

a_Ku(k-K)＝a_Kc+a_Ks(k-K)+a_Kv_n(k-K)

考虑到式(A7)，有如下

u(k)-a₁u(k-1)-...-a_Ku(k-K)＝c(1-a₁-..-a_K)+,k＝1,…M,M≥K+1 (A13)

(v_n(k)-a₁v_n(k-1)-..-a_Kv_n(k-K))

令

q＝c(1-a₁-..-a_K)

u＝[u(k),u(k+1),…u(k+M)]^T

p＝[q,a₁,…a_K]^T

推出方程和算法(4)～(6)。注意到(A13)式中的数量至少需要K+1个，这是因为需要被辨识的参数多于普通的AR过程。

与([4])中类似，如果传感器含有附加的高斯白噪声，估计器产生了的最优估计。

$E\left(e^2\right)=E((a_1-{\hat{a}}_1)u(k-1)+...+(a_K-{\hat{a}}_{K/2})u(k-K)-(q-\hat{q}){)}^2)---\left(A14\right)$

三、噪声协方差估计

众所周知，飞机以θ和γ进行坐标转换，转换矩阵从机体坐标系到地理坐标系如下给出：

$A_b^g=\left[\begin{array}{ccc}\cos \mathrm{\γ}\cos \mathrm{\φ}-\sin \mathrm{\γ}-\sin \mathrm{\φ}\sin \mathrm{\θ}& -\sin \mathrm{\φ}\cos \mathrm{\θ}& \cos \mathrm{\φ}\sin \mathrm{\γ}+\sin \mathrm{\φ}\sin \mathrm{\θ}\cos \mathrm{\γ}\\ \cos \mathrm{\γ}\sin \mathrm{\φ}+\sin \mathrm{\γ}\cos \mathrm{\φ}\sin \mathrm{\θ}& \cos \mathrm{\φ}\cos \mathrm{\θ}& \sin \mathrm{\γ}\sin \mathrm{\φ}-\cos \mathrm{\φ}\sin \mathrm{\θ}\cos \mathrm{\γ}\\ -\cos \mathrm{\θ}\sin \mathrm{\γ}& \sin \mathrm{\θ}& \cos \mathrm{\θ}\cos \mathrm{\γ}\end{array}\right]---\left(B1\right)$

磁场观测量转换到地理坐标系，如下表示

$B_g=A_b^g{B}_b---\left(B2\right)$

考虑姿态测量中的噪声，θ₀、γ₀表示当前值，Δθ和Δγ表示量测噪声。因此包含噪声的磁场可以如下表示

B_b＝B_b0+ΔB_b

因此地理坐标系中的磁场误差可以写成：

$\mathrm{\Δ}{B}_g\≈A_b^g\mathrm{\Δ}{B}_b+\mathrm{\Δ}{A}_b^g{B}_{b0}---\left(B3\right)$

公式中高于二阶的误差项被忽略不计，是姿态测量噪声导致的扰动，分析姿态扰动表达式可以获知姿态测量误差如何导致总场误差。总体而言，这些操作是繁琐和冗长的，忽略高于1阶的次数项，(3)右侧第二项可以被写成如下

$\mathrm{\Δ}{A}_b^g{B}_{n0}=C(B_{b0},{\mathrm{\φ}}_0,{\mathrm{\θ}}_0,{\mathrm{\λ}}_0)\left[\begin{array}{c}\mathrm{\Δ\φ}\\ \mathrm{\Δ\θ}\\ \mathrm{\Δ\γ}\end{array}\right]..---\left(B4\right)$

进一步地，载机在水平面以给定的航向、俯仰、横滚飞行，θ和γ都比较小,sinθ和sinγ可近似为θ和γ,并简化成如下：

从(3)和(4)看出总场误差B_b是高斯的。而在本文中传感器噪声和角度误差都认为是高斯型噪声。注意到传感器噪声独立于姿态误差，坐标转换后总误差的协方差矩阵可如下表示：

$\begin{array}{c}{R}_{B_g}=E\left\{\mathrm{\Δ}{B}_g\mathrm{\Δ}{B}_g^T\right\}=\mathrm{diag}[{\mathrm{\σ}}_{\mathrm{bx}}^2,{\mathrm{\σ}}_{\mathrm{by}}^2,{\mathrm{\σ}}_{\mathrm{bz}}^2]+\\ +E\left\{C(B_{g0},{\mathrm{\φ}}_0,{\mathrm{\θ}}_0,{\mathrm{\λ}}_0)\left[\begin{array}{c}\mathrm{\Δ\φ}\\ \mathrm{\Δ\θ}\\ \mathrm{\Δ\γ}\end{array}\right]\right[\mathrm{\Δ\φ},\mathrm{\Δ\θ},\mathrm{\Δ\γ}\left]C^T(B_{g0},{\mathrm{\φ}}_0,{\mathrm{\θ}}_0,{\mathrm{\λ}}_0)\right\}\end{array}---\left(B6\right)$

注意(6)式右侧第二项，它是独立于转换坐标系的。.

本发明提供了一个用于在中纬度地区进行地磁观测的空地一体的航空地磁场观测系统，由机载部分和移动地面站部分组成。首先，利用移动地面站采集的地磁场磁场信息，采用拟自回归(AR)过程建立局部地磁场模型，对该模型的参数进行识别；其次，在系统的机载任务处理单元中，利用地面识别的局部地磁场模型参数，构建了沿飞行航路的地磁场递推形式的动态模型；再次，根据机载部分的传感器配置和姿态传感器误差，构建了沿飞行航路的测量误差模型；最后，综合上述动态模型和测量误差模型，设计了一个在机载任务处理单元中执行的扩展Kalman滤波算法，提高了在地理坐标系中表示的地磁场向量的测量精度，其中的新息还可用于地磁场异常的检测。

1.本发明在局部地区，沿着与磁通量在水平面上投影共线的基线移动的基站，用于收集本地磁场强度向量和总量数据，利用所采集到的数据识别出一个适当阶的本地磁场拟AR模型的各个模态。

2.本发明所述扩展Kalman滤波器(EKF),使用机载光泵磁强计和三轴磁通门磁强计采集到的数据，利用1中识别出的当地磁场模型的参数，构造出沿预定航路上地磁场的动态模型，从姿态测量噪声和传感器噪声中，估计出在预先定义好的地理坐标系中表达的地磁场的状态。

3.本发明能够输出高精度的地磁场观测数据。

4.本发明通过一个磁异常检测器检测可能出现的磁异常。

尽管上文对本发明进行了详细说明，但是本发明不限于此，本技术领域技术人员可以根据本发明的原理进行各种修改。因此，凡按照本发明原理所作的修改，都应当理解为落入本发明的保护范围。

Claims (10)

1.一种空地一体的地磁场联合观测方法，其特征在于，包括：

在地磁场测量任务开展前，利用移动地面站采集的地磁场信息，确定局部地磁场模型的模型参数；

在执行地磁场测量任务时，利用机载部分在飞行航路上采集的地磁场信息和所述模型参数，得到沿飞行航路上的动态模型；

利用机载部分获取的载机按预定姿态飞行时的姿态测量误差，得到沿飞行航路上的测量误差模型；

利用所述动态模型和所述测量误差模型，得到飞行航路上采集的地磁场信息在地理坐标系中的估计。

2.根据权利要求1所述的方法，其特征在于，在待测区域内选取若干个观测点，利用移动地面站采集所述若干个观测点的地磁场信息，并利用所述若干个观测点的地磁场信息，估计所述待测区域内的局部地磁场模型的模型参数，其中，所述地磁场信息包括地磁场总量和向量，所述模型参数包括所述局部地磁场模型的直流分量和多个谐波分量的系数。

3.根据权利要求2所述的方法，其特征在于，所述若干个观测点均匀分布在一条基线上，所述基线是一条原点在参考点x₀，且与磁北方向重合的长度为L的直线。

4.根据权利要求1所述的方法，其特征在于，利用所述局部地磁场模型的模型参数，对所述机载部分在飞行航路上采集的地磁场信息进行分解，得到各个方向的直流分量和谐波分量，并将各个方向上的谐波分量进行离散化处理，得到所述飞行航路上的动态模型。

5.根据权利要求1所述方法，其特征在于，机载部分通过其定姿定位系统获取姿态测量误差。

6.根据权利要求1-5任意一项所述方法，其特征在于，利用所述动态模型和所述测量误差模型，得到在地理坐标系中的卡尔曼滤波器，并通过所述卡尔曼滤波器，得到飞行航路上采集的地磁场信息在地理坐标系中的估计。

7.根据权利要求6所述方法，其特征在于，还包括：利用所述卡尔曼滤波器内的数据流，检测地磁场是否发生异常。

8.一种空地一体的地磁场联合观测系统，其特征在于，包括：

移动地面站，用于在地磁场测量任务开展前，利用移动地面站采集的地磁场信息，确定局部地磁场模型的模型参数；

机载部分，用于在执行地磁场测量任务时，利用在飞行航路上采集的地磁场信息和所述模型参数，得到沿飞行航路上的动态模型，利用获取的载机按预定姿态飞行时的姿态测量误差，得到沿飞行航路上的测量误差模型，并利用所述动态模型和所述测量误差模型，得到飞行航路上采集的地磁场信息在地理坐标系中的估计。

9.根据权利要求8所述的系统，其特征在于，所述机载部分利用所述动态模型和所述测量误差模型，得到在地理坐标系中的卡尔曼滤波器，并通过所述卡尔曼滤波器，得到飞行航路上采集的地磁场信息在地理坐标系中的估计。

10.根据权利要求9所述的系统，其特征在于，还包括：

磁异常检测器，用于利用所述卡尔曼滤波器内的数据流，检测地磁场是否发生异常。