一种时间域航空电磁探测系统接收线圈运动噪声的数值仿真方法与流程

本发明涉及一种时间域航空电磁探测的资料处理方法，尤其是接收线圈运动噪声的数值仿真方法。

背景技术：

航空电磁探测作为一种高效快速的物探方法，是利用飞机承载探测仪器，在航行过程中观测地球物理场发生的各种变化，并根据电磁感应原理获取大地电阻率信息，用于勘查地下矿产等资源或进行地质问题的研究。在时间域航空电磁数据处理中，噪声是影响航空电磁探测深度与异常分辨能力的重要因素。其中，航空电磁探测系统在飞行测量过程中，由于飞行环境、飞机姿态、飞行速度等原因，引起接收线圈在地磁场中摆动，切割地磁场磁力线进而产生感应电动势，即为线圈运动噪声。因其具有频率低、幅度大的特点而成为航空电磁探测中的主要噪声源之一。线圈运动噪声会导致观测数据与系统平稳飞行结果之间存在较大偏差，因而影响了观测数据的一致性，对航空电磁探测成像结果造成干扰。

munkholm(motion-inducednoisefromvibrationofamovingtemdetectorcoil:characterizationandsuppression.journalofappliedgeophysics,1997)利用三分量磁场之和在主场方向上的投影与运动噪声之间耦合最小的原理抑制运动噪声；

davis(pendulummotioninairbornehemsystems.explorationgeophysics,2006)通过分析运动噪声振荡周期特点设计滤波器进行压制等；

吉林大学学报(地球科学版),2013(时间域航空电磁数据线圈运动噪声去除方法仿真研究),采用多项式拟合全波半周期运动噪声，运用拉格朗日优化算法求解运动噪声多项式系数，然后在半周期数据中去除运动噪声。

cn104865608a公开了一种“时间域航空电磁法运动噪声检测装置及抑制方法”。该方法设计一种运动噪声检测装置，包括地面和空中两部分，空中由直升机或飞艇搭载，地面是装有三个磁通门传感器的箱体和地面数据收录设备，根据采集数据对运动噪声进行检测并加以抑制。上述关于时间域航空电磁探测系统接收线圈运动噪声的压制方案在一定程度上降低了运动噪声的干扰，但是仍然缺乏较为完善的针对运动噪声的原理及计算方法的研究。到目前为止，尚未见基于线圈姿态的时间域航空电磁探测系统接收线圈运动噪声的计算方法。

技术实现要素：

本发明的目的就在于针对上述现有技术的不足，提出一种基于线圈姿态的时间域航空电磁探测系统接收线圈运动噪声的计算方法。该方法针对时间域航空电磁探测系统接收线圈摆动形成的噪声与干扰问题，基于系统的线圈姿态测量，首先定义了接收线圈运动噪声的概念，依据飞行测量过程中接收线圈切割地磁场磁力线产生感应电动势推导得出运动噪声计算公式，其次利用飞行测量过程中的地理信息选取地磁场参数，并建立相对应的惯性坐标系和接收线圈坐标系，基于线圈姿态建立旋转矩阵以进行坐标转换，最终利用运动噪声计算公式得出数值仿真结果。

一种时间域航空电磁探测系统接收线圈运动噪声的数值仿真方法，包括以下步骤：

a、录入航空电磁辅助信息数据；

b、依据航空电磁探测地点及时间选取地磁场参数；

c、根据法拉第电磁感应定律推导航空电磁探测接收线圈运动噪声理论计算公式；

d、分别建立惯性坐标系和接收线圈坐标系并定义法向量；

e、依据惯性坐标系与接收线圈坐标系关系，引入接收线圈姿态的旋转矩阵d；

f、利用旋转矩阵d推导出非平稳状态下接收线圈运动噪声计算公式；

g、根据辅助信息数据及地磁场参数计算各分量线圈运动噪声εi；

h、判断是否录入计算完成，否，返回到上一步；

i、是，输出结果。

步骤a所述的航空电磁辅助信息数据，包括飞行测量过程中接收线圈的摆动角度信息和飞行方向信息。

步骤b所述的地磁场参数包括：磁场总分量f、垂直分量z和北向分量x。

首先，录入航空电磁辅助信息数据，得到接收线圈各姿态角度数据，分别是摇摆角度数据α、俯仰角度数据β和偏航角度数据γ，并确定航空电磁探测飞行时间、飞行地点及航向θ。

其次，依据航空电磁探测飞行时间及地点，在国家地磁台网中心下载相应地磁场参数，分别是磁场总分量f、垂直分量z和北向分量x。

根据得到的地磁场参数，进一步计算东向分量y：

y＝-h·sin(d),(1)

其中，为地磁场磁偏角，h＝f·cos(i)为地磁场水平分量，而地磁场磁倾角i通过关系式：计算。

再依据航向将计算所得的原始地磁场分量表示为航空电磁探测系统飞行方向上的地磁场分量be＝[bxbybz]，计算公式如下：

bx＝x·cosθ+y·sinθ,

by＝-x·sinθ+y·cosθ,(2)

bz＝z.

第三步，根据法拉第电磁感应定律，分析线圈运动噪声成因，推导出航空电磁探测接收线圈运动噪声理论计算公式。

运动噪声是在航空电磁探测过程中接收线圈切割地磁场磁力线，引起接收线圈内磁通量变化而产生的感应电动势。根据法拉第电磁感应定律，在飞行测量中接收线圈所产生的感应电动势可表示为：

其中，φ为磁通量。通过分析在测线上飞行探测时感应的场源类型，可以得出场源主要包括地磁场be和二次场ba。因此公式(3)可以表示为：

其中，s为接收线圈面积；地磁场be为常数，其导函数为零，因此公式(4)中的第二项为零；二次场ba的数量级远小于地磁场be，故公式(4)中的第三项可以忽略；此外，公式(4)中第四项表示接收线圈接收到的二次场响应，为航空电磁探测中的待测信号。因此，考虑到接收线圈内感应电动势的理论分析以及运动噪声的成因，可以定义运动噪声的理论计算公式为：

第四步，建立接收线圈坐标系，该坐标系用于表示飞行过程中接收线圈摆动后的坐标(x′,y′,z′)。

设系统在不平稳状态下飞行，x′轴平行飞行方向，z′轴垂直于接收线圈平面。接收线圈坐标系为笛卡尔坐标系，即线圈旋转满足右手螺旋定则。定义a′为接收线圈平面的法向量，用于表示接收线圈姿态与线圈面积的关系。设定法向量a为(1,0,0)^t,(0,1,0)^t,(0,0,1)^t来分别表示平稳状态下x、y、z方向的接收线圈法向量。

第五步，依据惯性坐标系与接收线圈坐标系关系，引入基于接收线圈姿态的旋转矩阵d。

定义接收线圈旋转方向为摇摆、俯仰、偏航。惯性坐标系与接收线圈坐标系关系可表示为下式：

a＝dγdβdαvb＝da′,(6)

其中，a为惯性坐标系法向量；dα、dβ、dγ分别表示接收线圈坐标系下的摇摆、俯仰和偏航姿态的旋转矩阵，d则为总旋转矩阵。各旋转矩阵表达式如下：

第六步，推导非平稳状态下接收线圈运动噪声计算公式。

在非平稳状态下接收线圈面积s是矢量，且表示为接收线圈的有效面积与非平稳状态下接收线圈法向量的乘积，即s＝sa′。同时，非平稳状态下接收线圈法向量a′满足：a′＝d′a，其中，d′为旋转矩阵d的逆矩阵。因此，上述公式带入公式(5)中，可以得到：

公式(9)即为航空电磁探测接收线圈运动噪声理论计算方法。

第七步，带入第一步导入的辅助信息数据及第二步分解的地磁场参数计算各分量接收线圈运动噪声εi(i＝x,y,z)，根据带入的数据长度进行判断，直至计算完成，输出该分量接收线圈运动噪声计算结果。

有益效果：本发明从理论上有效地数值仿真了接收线圈在地磁场中摆动而引入的运动噪声，对计算的接收线圈运动噪声进行航空电磁数据预处理，并与实测数据进行对比分析。经分析验证，计算结果既真实又可靠。本方法直接利用接收线圈辅助信息和地磁场参数进行计算，原理清晰，易于实施。经试验，本方法对时间域航空电磁探测接收线圈运动噪声的计算精度较高，为进一步去除运动噪声对航空电磁数据的干扰提供了技术支持。按照本发明的方法，无论接收线圈摆动角度为多少，均可进行运动噪声计算，且计算结果精确。

附图说明

图1是时间域航空电磁探测系统接收线圈运动噪声数值仿真方法流程图。

图2是录入的航空电磁辅助信息数据中的姿态角度数据。

图中实线为roll姿态角度数据，---虚线为pitch姿态角度数据，……线为yaw姿态角度数据。

图3是数值仿真得到的z分量接收线圈运动噪声。

具体实施方案

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

一种时间域航空电磁探测系统接收线圈运动噪声的数值仿真方法，包括以下步骤：

a、录入航空电磁辅助信息数据；

b、依据航空电磁探测地点及时间选取地磁场参数；

c、根据法拉第电磁感应定律推导航空电磁探测接收线圈运动噪声理论计算公式；

d、分别建立惯性坐标系和接收线圈坐标系并定义法向量；

e、依据惯性坐标系与接收线圈坐标系关系，引入接收线圈姿态的旋转矩阵d；

f、利用旋转矩阵d推导出非平稳状态下接收线圈运动噪声计算公式；

g、根据辅助信息数据及地磁场参数计算各分量线圈运动噪声εi；

h、判断是否录入计算完成，否，返回到上一步；

i、是，输出结果。

以z分量线圈运动噪声计算为例，计算数据采用项目组于2015年9月18日在湖南省岳阳临湘市桃矿进行野外飞行实验中采集到的航空电磁辅助信息数据。

首先，录入航空电磁辅助信息数据，得到接收线圈各姿态角度数据，分别是摇摆(roll)角度数据α、俯仰(pitch)角度数据β和偏航(yaw)角度数据γ，如图2所示，并确定航空电磁探测飞行时间、飞行地点及航向θ＝5°(北偏西)。

其次，依据航空电磁探测飞行时间及地点，在国家地磁台网中心下载2015年9月湖南省地磁场参数，分别是磁场总分量f＝49549nt、垂直分量z＝36092nt和北向分量x＝33855nt。

根据得到的地磁场参数，进一步计算东向分量y：

y＝-h·sin(d)＝-2512nt,(1)

其中，地磁场磁偏角地磁场水平分量h＝f·cos(i)＝33948nt，地磁场磁倾角

再依据航向将计算所得的原始地磁场分量表示为航空电磁探测系统飞行方向上的地磁场分量，计算公式如下：

bx＝x·cosθ+y·sinθ,

by＝-x·sinθ+y·cosθ,(2)

bz＝z.

得到系统飞行方向上的地磁场分量

be＝[bxbybz]＝[33507-5452136092]。

第三步，根据法拉第电磁感应定律，分析线圈运动噪声成因，推导出航空电磁探测接收线圈运动噪声理论计算公式。

运动噪声是在航空电磁探测过程中接收线圈切割地磁场磁力线，引起接收线圈内磁通量变化而产生的感应电动势。根据法拉第电磁感应定律，在飞行测量中接收线圈所产生的感应电动势可表示为：

其中，φ为磁通量。通过分析在测线上飞行探测时感应的场源类型，可以得出场源主要包括地磁场be和二次场ba。因此公式(3)可以表示为：

其中，s为接收线圈面积；地磁场be为常数，其导函数为零，因此公式(4)中的第二项为零；二次场ba的数量级远小于地磁场be，故公式(4)中的第三项可以忽略；此外，公式(4)中第四项表示接收线圈接收到的二次场响应，为航空电磁探测中的待测信号。因此，考虑到接收线圈内感应电动势的理论分析以及运动噪声的成因，可以定义运动噪声的理论计算公式为：

第四步，建立接收线圈坐标系，该坐标系用于表示飞行过程中接收线圈摆动后的坐标(x′,y′,z′)。

设系统在不平稳状态下飞行，x′轴平行飞行方向，z′轴垂直于接收线圈平面。接收线圈坐标系为笛卡尔坐标系，即线圈旋转满足右手螺旋定则。定义a′为接收线圈平面的法向量，用于表示接收线圈姿态与线圈面积的关系。设定法向量a为(0,0,1)^t来表示平稳状态下z方向的接收线圈法向量。

第五步，依据惯性坐标系与接收线圈坐标系关系，引入基于接收线圈姿态的旋转矩阵d。

定义接收线圈旋转方向为摇摆、俯仰、偏航。惯性坐标系与接收线圈坐标系关系可表示为下式：

a＝dγdβdαvb＝da′,(6)

其中，a为惯性坐标系法向量；dα、dβ、dγ分别表示接收线圈坐标系下的摇摆、俯仰和偏航姿态的旋转矩阵，d则为总旋转矩阵。各旋转矩阵表达式如下：

在非平稳状态下接收线圈面积s是矢量，且表示为接收线圈的有效面积与非平稳状态下接收线圈法向量的乘积，即s＝sa′(接收线圈面积单位化，即s＝1)。同时，非平稳状态下接收线圈法向量a′满足：a′＝d′a，其中，d′为旋转矩阵d的逆矩阵。因此，上述公式带入公式(5)中，可以得到：

公式(9)即为z分量航空电磁探测接收线圈运动噪声理论计算方法。

第七步，带入第一步导入的辅助信息数据及第二步分解的地磁场参数计算z分量线圈运动噪声εz，根据带入的数据长度进行判断，直至计算完成，输出z分量接收线圈运动噪声计算结果如图3所示。