CN108984818A - 固定翼时间域航空电磁数据拟三维空间约束整体反演方法 - Google Patents

固定翼时间域航空电磁数据拟三维空间约束整体反演方法 Download PDF

Info

Publication number
CN108984818A
CN108984818A CN201810492825.8A CN201810492825A CN108984818A CN 108984818 A CN108984818 A CN 108984818A CN 201810492825 A CN201810492825 A CN 201810492825A CN 108984818 A CN108984818 A CN 108984818A
Authority
CN
China
Prior art keywords
inversion
attitude angle
dimensional space
model
conductivity
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN201810492825.8A
Other languages
English (en)
Inventor
朱凯光
张琼
彭聪
王昊
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Jilin University
Original Assignee
Jilin University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Jilin University filed Critical Jilin University
Priority to CN201810492825.8A priority Critical patent/CN108984818A/zh
Publication of CN108984818A publication Critical patent/CN108984818A/zh
Pending legal-status Critical Current

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V3/00—Electric or magnetic prospecting or detecting; Measuring magnetic field characteristics of the earth, e.g. declination, deviation
    • G01V3/08—Electric or magnetic prospecting or detecting; Measuring magnetic field characteristics of the earth, e.g. declination, deviation operating with magnetic or electric fields produced or modified by objects or geological structures or by detecting devices
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V3/00—Electric or magnetic prospecting or detecting; Measuring magnetic field characteristics of the earth, e.g. declination, deviation
    • G01V3/38—Processing data, e.g. for analysis, for interpretation, for correction
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/15—Correlation function computation including computation of convolution operations
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/15—Vehicle, aircraft or watercraft design
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Computational Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • Remote Sensing (AREA)
  • Geology (AREA)
  • Computing Systems (AREA)
  • Algebra (AREA)
  • Geophysics (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • Environmental & Geological Engineering (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Electromagnetism (AREA)
  • Automation & Control Theory (AREA)
  • Aviation & Aerospace Engineering (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本发明涉及一种固定翼时间域航空电磁数据拟三维空间约束整体反演方法，首先设置初始模型和正则化系数，对各测点数据进行正演计算测区数据误差项；其次，若误差项未达到反演终止条件，建立地下各层电导率和姿态角度的拟三维空间约束整体反演目标函数；最后再次判断数据误差项是否达到反演终止条件，即为反演结果。本发明不仅考虑了各测线间的连续性，而且加入了飞行姿态角度。引入拟三维空间约束反演确保反演结果精确；不仅对单测点地下各层电阻率和单测线各测点反演参数进行约束，且对测区各测线进行整体约束；避免了因姿态角度不精确或者缺失无法直接进行姿态角度校正。拟三维空间约束整体反演结果优于拟二维整体反演结果。

Description

固定翼时间域航空电磁数据拟三维空间约束整体反演方法
技术领域
本发明涉及一种时间域航空地球物理勘探方法，尤其是固定翼时间域航空电磁数据拟三维空间约束整体反演方法。
背景技术
固定翼航空电磁探测具有工作效率高、探测面积广以及探测深度深等优点，在矿产探测、地质普查等方面具有广泛的应用。同时，固定翼航空电磁探测存在数据量大、测区数据飞行姿态不一致等问题，导致数据解释困难，反演结果成像不完整。
航空电磁一维反演对单测点数据进行独立反演，忽略了相邻测点间的相关性，导致反演结果不连续，测线或测区整体反演结果存在较大的误差。航空电磁数据二维反演结果和三维反演结果更能准确描述地下地质分布情况，但反演算法复杂，计算量大，不仅对计算机内存和速度要求较高，且耗时过长，在野外实测数据中难以实施。综合航空电磁一维反演和二维反演、三维反演的优势，国内外研究了基于一维反演的约束反演。横向约束反演基于测线上邻近测点数据的连续性，对相邻测点反演参数进行相互约束，保证了反演结果沿剖面方向的平滑性和连续性。
毛立峰等公开了《飞行高度同时反演的固定翼航空瞬变电磁一维反演》.以层状模型的固定翼时间域航空电磁多分量理论响应数据为例，将飞行高度作为一个待反演参数与电阻率参数一并反演，但并未考虑各飞行姿态角度的影响。到目前为止，仍未见利用footprint区域内测点进行约束且把飞行姿态角度作为反演参数的拟三维空间约束整体反演。地球物理学报,2011,8,P2136-2147.
殷长春等公开了航空电磁单个测点“敏感区域”通常为飞行高度的3-4倍，仅对邻近测点进行约束可能造成有效信息的丢失，造成测区各测线结果不一致，导致整体模型不够平滑，但并未见在拟三维反演中根据footprint进行相互约束。《航空电磁勘查技术发展现状及展望》.地球物理学报,2015,58(8):P2637-2653.
王琦等公开了《固定翼时间域航空电磁数据拟二维整体反演》.基于一维Tikhonov正则化反演理论，通过对测线上相邻测点进行约束，实现了对整条测线数据同时反演的固定翼航空电磁数据的拟二维整体反演算法，但并未对测线间数据进行约束，测区反演结果会出现测线间不连续的现象。地球物理学进展,2016,3:P1173-1180.
此外，在实际探测中，固定翼航空电磁系统不能持续保持在平稳状态下飞行，根据已知资料，接收线圈的pitch姿态角度在10s内可以从-20°变为20°，这将会使数据变化30.47％至51.3％。当线圈产生姿态角度时，不能采用平稳飞行时的正演公式，若忽略姿态角度的影响，将严重影响反演结果。
CN106338774A公开了“一种基于电导率-深度成像的时间域航空电磁数据反演方法”。该方法首先采用基于查表法的电导率-深度成像方法获得地下介质的视电导率和视深度，然后以此构建反演初始模型，最后应用阻尼特征参数法进行反演完成时间域航空电磁数据的组合解释。该方法在获得近似成像结果的基础上进行反演，但是并未涉及footprint确定反演约束范围和飞行姿态角度的反演。
发明内容
本发明的目的在于针对上述现有技术的不足，提供一种准确的固定翼时间域航空电磁数据拟三维空间约束整体反演方法。
本发明主要是针对现有反演方法仅对邻近测点进行约束，忽略了航空电磁相邻测线数据的连续性，可能造成有效信息的丢失；此外现有反演方法忽略了飞行姿态角度(发现发射线圈、接收线圈的俯仰姿态和吊舱同向摆动)对反演结果的影响。
本发明的思想是：首先设置初始模型和正则化系数，对各测点数据进行正演计算，并利用正演结果计算测区数据误差项；其次，若误差项未达到反演终止条件，建立关于地下各层电导率和姿态角度的拟三维空间约束整体反演目标函数，求解迭代方程，得到新反演模型；最后再次判断数据误差项是否达到反演终止条件，若达到反演终止条件，此反演模型即为反演结果，否则，继续迭代求解。
固定翼时间域航空电磁数据拟三维空间约束整体反演方法，包括以下步骤：
a、录入固定翼时间域航空电磁数据对应的先验信息，本文采用电导率深度成像结果作为模型先验信息，包括地下各层电导率和飞行姿态角度；
b、设置待反演参数包括地下各层电导率和飞行姿态角度(发射线圈俯仰姿态角度、接收线圈俯仰姿态角度和同向摆动角度)初始模型，为加速收敛过程，本文把电导率深度成像结果设为电导率初始模型，包括地下各层电阻率模型和飞行姿态角度；
c、对各测点数据进行考虑发射线圈俯仰姿态、接收线圈俯仰姿态和吊舱同向摆动的正演计算，在飞机飞行过程中，吊舱的摆动会影响系统的几何参数，发射线圈旋转存在姿态角度时，发射线圈的发射磁矩会发生改变，接收线圈旋转存在姿态角度时，接收分量的方向会发生改变；为得到时域电磁响应，首先将频域二次场变换到s域，再通过拉普拉斯逆变换，得到时域感应电动势；
d、利用正演结果计算数据误差项，判断是否达到反演终止条件，是，输出反演结果；否，进行下一步；
e、建立关于测区测点电导率和飞行姿态角度(发射线圈俯仰姿态角度、接收线圈俯仰姿态角度和同向摆动角度)的拟三维空间约束整体反演目标函数，将发射线圈的俯仰姿态角度、接收线圈的俯仰姿态角度和吊舱的同向摆动角度引入拟三维空间约束整体反演中，此时待求参数包含地下各层电导率和系统飞行姿态角度；
f、利用超松弛迭代预处理共轭梯度算法求解迭代方程，得到新反演模型，令目标函数对待求反演模型参数的偏导等于零，利用超松弛迭代预处理共轭梯度算法求解目标函数的最优解；
g、利用新反演模型计算数据误差项，并判断是否达到反演终止条件，否，返回到步骤e，建立关于测区测点电导率和姿态角度的新的拟三维空间约束整体反演目标函数；
h、是，输出反演结果。
有益效果：本发明公开的固定翼时间域航空电磁数据拟三维空间约束整体反演方法，(1)将飞行姿态角度引入拟三维空间约束反演中，可以保证反演结果精确性；(2)本发明不仅对单测点地下各层电阻率和单测线各测点反演参数进行约束，且对测区各测线进行约束，使测区反演结果平滑；(3)本发明避免了姿态角度不精确或者缺失，无法直接进行姿态角度校正，造成的反演结果的误差。
本发明不仅考虑了各测线间的连续性，而且加入了飞行姿态角度，使整体反演过程更加的精确。
从反演结果对比可以看出，拟三维空间约束整体反演结果优于拟二维整体反演结果，尤其在存在噪声的情况下。
附图说明
图1固定翼时间域航空电磁数据拟三维空间约束整体反演方法流程图。
图2是大地模型示意图。
图3是拟二维整体反演结果
图3a第一条测线反演结果；
图3b第二条测线反演结果；
图3c第三条测线反演结果；
图3d第四条测线反演结果；
图3e第五条测线反演结果；
图3f第六条测线反演结果。
图4是拟三维空间约束整体反演结果。
图4a第一条测线反演结果；
图4b第二条测线反演结果；
图4c第三条测线反演结果；
图4d第四条测线反演结果；
图4e第五条测线反演结果；
图4f第六条测线反演结果。
具体实施方式
下面结合附图和实施例对本发明作进一步的详细说明：
固定翼时间域航空电磁数据拟三维空间约束整体反演方法，包括以下步骤：
a、录入固定翼时间域航空电磁数据对应的先验信息，
所述的录入固定翼时间域航空电磁数据对应的先验信息，包括地下各层电导率和飞行姿态角度，整理数据为RD(L×N×(M+3),1)，其中M为反演大地层数；
所述的录入固定翼时间域航空电磁数据，设待反演时间域航空电磁测区数据D包含L条测线，每条测线有N个测点，每个测点有t个时间道，整理数据为D(L×N×t,1)，
其中，T为转置符号，测点数据
b、设置待反演参数包括地下各层电导率和飞行姿态角度(发射线圈俯仰姿态角度、接收线圈俯仰姿态角度和同向摆动角度)初始模型，为加速收敛过程，本文把电导率深度成像结果设为电导率初始模型，包括地下各层电阻率模型和飞行姿态角度；
所述的设置待反演参数包括地下各层电导率和飞行姿态角度(发射线圈俯仰姿态角度、接收线圈俯仰姿态角度和同向摆动角度)初始模型，为加速收敛过程，本文把电导率深度成像结果设为电导率初始模型，包括地下各层电阻率模型和飞行姿态角度。
c、对各测点数据进行考虑发射线圈俯仰姿态、接收线圈俯仰姿态和吊舱同向摆动的正演计算，在飞机飞行过程中，吊舱的摆动会影响系统的几何参数，发射线圈旋转存在姿态角度时，发射线圈的发射磁矩会发生改变，接收线圈旋转存在姿态角度时，接收分量的方向会发生改变；为得到时域电磁响应，首先将频域二次场变换到s域，再通过拉普拉斯逆变换，得到时域感应电动势；
所述的对各测点数据进行考虑发射线圈俯仰姿态、接收线圈俯仰姿态和吊舱同向摆动的正演计算，飞行方向设为x轴方向，发射线圈中心在地面的投影设为坐标系的原点，发射线圈中心坐标为(0,0,ht)，接收线圈中心坐标为(x0,y0,z0)，系统水平收发距为r＝√(x0 2+y0 2)。设地下为N层各向同性均匀的大地模型，层状大地的第N层电导率和厚度分别为σN(S/m)、dN(m)。
在飞机飞行过程中，吊舱的摆动会影响系统的几何参数，设吊舱同向摆动角度为γ，接收线圈坐标为其中，γ0为平稳飞行下吊舱悬绳与z轴负方向的夹角，此时，系统水平收发距为r'＝√(x'0 2+y'0 2)。其中，AR为接收线圈有效面积，当发射线圈内通以角频率为ω，强度为ejωt的电流时，接收线圈(x'0,y'0,z'0)处的二次场频域表达式为，
H(x'0,y'0,z'0,ω)＝G·m， (2)
其中，为发射磁矩，m0为磁矩系数，为单位方向矢量，若仅在z方向存在发射，发射线圈旋转存在姿态角度时，发射线圈的发射磁矩会发生改变，设发射线圈的俯仰姿态角度为α，产生旋转矩阵系统实际发射磁矩m'＝RT·m；二次场格林张量为
λ为积分变量，反射系数R0通过递归方法计算，J0和J1为零阶和一阶贝塞尔函数，含有贝赛尔函数的无限积分为汉克尔变换，可根据Guptasarma滤波计算。接收线圈旋转存在姿态角度时，接收分量的方向会发生改变，设接收线圈的俯仰姿态角度为β，产生旋转矩阵系统实际接收到二次场为H'＝RR T·G·RT·m。
设发射电流为阶跃电流，为得到时域电磁响应，首先将频域二次场变换到s域，再通过拉普拉斯逆变换，得到时域感应电动势为，
V(t)＝L-1[V(s)]＝ARμ0L-1[H(s)]＝ARμ0L-1[H'(ω)|jω＝s] (3)
μ0＝4π×10-7H/m为空气磁导率，式(3)中的拉普拉斯逆变换通过G-S变换计算。
d、利用正演结果计算数据误差项，判断是否达到反演终止条件，是，输出反演结果；否，进行下一步；
所述的利用正演结果计算数据误差项Φd(m)，
其中D为测区航空电磁数据，如步骤a所述；fforward(m)为模型的正演响应；为测区航空电磁数据估计的误差协方差矩阵。
判断是否达到反演终止条件，反演的终止条件为，
Φd(m)/Nd≤1 (5)
其中，Nd为反演参数总个数，Nd＝(M+3)*N*L，
如果此时的反演参数m满足式(5)，反演结束，输出反演结果；否，进入下一步。
e、建立关于测区测点电导率和飞行姿态角度(发射线圈俯仰姿态角度、接收线圈俯仰姿态角度和同向摆动角度)的拟三维空间约束整体反演目标函数，将发射线圈的俯仰姿态角度、接收线圈的俯仰姿态角度和吊舱的同向摆动
角度引入拟三维空间约束整体反演中，此时待求参数包含地下各层电导率和系统飞行姿态角度；
所述的建立关于测区测点电导率和飞行姿态角度(发射线圈俯仰姿态角度、接收线圈俯仰姿态角度和同向摆动角度)的拟三维空间约束整体反演目标函数，将发射线圈的俯仰姿态角度、接收线圈的俯仰姿态角度和吊舱的同向摆动角度引入拟三维空间约束整体反演中，此时待求参数包含地下各层电导率和系统飞行姿态角度，
其中，为第L条测线第N个测点第M层的电导率，为第L条测线第N个测点的发射线圈俯仰姿态角度，为第L条测线第N个测点的接收线圈俯仰姿态角度，为第L条测线第N个测点的吊舱同向摆动角度。
空间约束整体反演算法的目标函数Φ(m)由数据误差项Φd(m)，模型先验信息误差项Φr(m)，模型纵向粗糙度Φv(m)，模型横向粗糙度Φh(m)和模型空间粗糙度Φl(m)构成：
Φ(m)＝Φd(m)+λ[λrΦr(m)+λvΦv(mσ)+λhΦh(m)+λlΦl(m)] (7)
λ为正则化因子，在反演迭代过程中通过线性搜索自适应迭代自动确定；λr、λh、λv和λl分别为模型先验信息、模型横向粗糙度、模型纵向粗糙度和模型空间粗糙度的正则化系数，根据地质情况，在初始反演时由人工决定，迭代过程中并不改变。
模型先验信息误差项mr为测区航空电磁数据对应的模型先验信息，如步骤a所述，本文采用电导率深度成像结果作为模型先验信息，为先验信息不确定度的方差矩阵。
模型纵向粗糙度由于模型纵向粗糙度仅对单个测点的各层电导率之间进行约束，此时只取m中电导率部分记为mσ，
Qv为模型纵向粗糙度矩阵，对单测点的各层电导率进行约束，本文根据该测点相邻层电导率的二阶差分计算，即
其中，为第L条测线第N个测点第m层的厚度，可得模型纵向粗糙度矩阵为：
模型横向粗糙度采用二阶导数计算模型参数的横向粗糙度，对相邻测点间的电导率和飞行姿态角度进行约束，则横向粗糙度矩阵为，
其中，Δx为相邻测点间的距离。模型横向粗糙度保证了单测线上相邻测点间的连续性，使反演结果沿测线平滑。
模型空间粗糙度同样采用二阶导数计算模型参数的空间粗糙度，对相邻测线间的电导率和飞行姿态角度进行约束，则横向粗糙度矩阵为，
其中，Δy为相邻测线间的距离。模型空间粗糙度确保了相邻测线间的连续性和反演结果测线间的平滑性。
f、利用超松弛迭代预处理共轭梯度算法求解迭代方程，得到新反演模型，令目标函数对待求反演模型参数的偏导等于零，利用超松弛迭代预处理共轭梯度算法求解目标函数的最优解；
所述利用超松弛迭代预处理共轭梯度算法求解迭代方程，得到新反演模型，令目标函数对待求反演模型参数的偏导等于零，求解目标函数的最优解，在第n次迭代过程为，
其中，数据误差项对待求反演模型参数的偏导为
Gn为雅克比矩阵，为第i个测点对应反演模型的正演数据fforward(mn)对第j个反演模型参数的偏数；模型先验信息误差项对待求反演模型参数的偏导为模型纵向粗糙度对待求反演模型参数的偏导为模型横向粗糙度对待求反演模型参数的偏导为模型空间粗糙度对待求反演模型参数的偏导为整合式(10)为Amn＝b的形式，
利用超松弛迭代预处理共轭梯度算法求解式(13)线性方程组，得到第n次迭代的反演模型参数mn。
g、利用新反演模型计算数据误差项，并判断是否达到反演终止条件，否，返回到步骤e，建立关于测区测点电导率和姿态角度的新的拟三维空间约束整体反演目标函数；
所述的利用新反演模型计算数据误差项，并判断是否达到反演终止条件，利用式(4)计算步骤f得到的反演模型参数mn对应的数据误差项，并带入式(5)，判断是否达到反演终止条件。如满足反演终止条件，进入步骤h，否，回到步骤e，利用步骤f得到的反演模型参数mn建立关于测区测点电导率和姿态角度的新的拟三维空间约束整体反演目标函数。
h、是，输出反演结果。
实施例
以一组典型的仿真数据为例，对拟三维空间约束整体反演的结果进行分析。
固定翼系统参数设置为：z分量发射线圈距离地面高度h＝120m，接收线圈位于(-130,0,70)m，发射线圈和接收线圈水平收发距r＝130m，对发射线圈面积、接收线圈面积、发射电流强度进行归一化，即ATX＝1m2，ARX＝1m2，阶跃发射电流强度I＝1A，每个测点off-time数据道时间为0.01ms～10ms之间等对数间隔的14个采样时刻。
大地模型参数设置为：整个测区6条测线，测线间距200m，测线总长500m，相邻测点间距10m。在第2条和第5条测线的101m-400m处存在一个地下埋深200m厚度为40m电导率为0.1S/m的异常体；在第3条和第4条测线的101m-150m和351m-400m处存在一个地下埋深200m厚度为40m电导率为0.1S/m的异常体，同时在151m-350m处存在一个地下埋深60m厚度为100m电导率为0.1S/m的异常体和一个地下埋深200m厚度为40m电导率为0.1S/m的异常体，半空间模型的电导率为0.01S/m，如图2所示。
为仿真实测数据，在电磁响应中按式在大地模型中加入高斯白噪声，其中，Gaussian(0,1)为方差为1均值为0的高斯分布，b为1ms时航空电磁响应的噪声，t为各时间道对应的抽道时间。
步骤一录入固定翼时间域航空电磁数据，待反演时间域航空电磁测区数据D包含6条测线，每条测线有50个测点，每个测点有14个时间道，整理数据为D(6×50×14,1)，
测点数据
录入固定翼时间域航空电磁数据对应的先验信息，拟三维整体反演采用15层定厚度大地模型，各层厚度均设为20m，本文采用电导率深度成像结果作为先验信息，包括地下15层各层电阻率模型和飞行姿态角度(发射、接收线圈俯仰姿态角度和同向摆动角度)，同样整理数据为RD(6×50×(15+3),1)。
步骤二设置待反演参数地下各层电导率和飞行姿态角度(发射线圈俯仰姿态角度、接收线圈俯仰姿态角度和同向摆动角度)初始模型，包括各测点地下15层各层电阻率模型和飞行姿态角度。
步骤三对各测点数据进行考虑发射线圈俯仰姿态、接收线圈俯仰姿态和吊舱同向摆动的正演计算。飞行方向设为x轴方向，发射线圈中心在地面的投影设为坐标系的原点，发射线圈中心坐标为(0,0,120)，接收线圈中心坐标为(-130,0,70)，系统水平收发距为r＝130。平稳飞行下吊舱悬绳与z轴负方向的夹角吊舱同向摆动角度为γ，接收线圈坐标为此时，系统水平收发距为r'＝√(x'0 2+y'0 2)。当发射线圈内通以角频率为ω，强度为ejωt的电流时，接收线圈(x'0,y'0,z'0)处的二次场频域表达式为，
H(x'0,y'0,z'0,ω)＝G·m， (2)
其中，为发射磁矩，m0为磁矩系数，仅在z方向存在发射，发射线圈旋转存在姿态角度时，发射线圈的发射磁矩会发生改变，设发射线圈的俯仰姿态角度为α，产生旋转矩阵系统实际发射磁矩m'＝RT·m；二次场格林张量为
λ为积分变量，反射系数R0通过递归方法计算，J0和J1为零阶和一阶贝塞尔函数，含有贝赛尔函数的无限积分为汉克尔变换，可根据Guptasarma滤波计算。接收线圈旋转存在姿态角度时，接收分量的方向会发生改变，设接收线圈的俯仰姿态角度为β，产生旋转矩阵系统实际接收到二次场为H'＝RRT·G·RT·m。
设发射电流为阶跃电流，为得到时域电磁响应，首先将频域二次场变换到s域，再通过拉普拉斯逆变换，得到时域感应电动势为，
V(t)＝L-1[V(s)]＝4π×10-7×L-1[H(s)]＝4π×10-7×L-1[H'(ω)|jω＝s] (3)
式(3)中的拉普拉斯逆变换通过G-S变换计算(Knight,1982)。
步骤四利用正演结果计算数据误差项Φd(m)，
其中D为测区航空电磁数据，如步骤一所述；fforward(m)为模型的正演响应；为测区航空电磁数据估计的误差协方差矩阵。
判断是否达到反演终止条件，反演的终止条件为，
Φd(m)/((15+3)*50*6)≤1 (5)
此时的反演参数m不满足式(5)，否，进入下一步。
步骤五建立关于测区测点电导率和飞行姿态角度(发射线圈俯仰姿态角度、接收线圈俯仰姿态角度和同向摆动角度)的拟三维空间约束整体反演目标函数，将发射线圈的俯仰姿态角度、接收线圈的俯仰姿态角度和吊舱的同向摆动角度引入拟三维空间约束整体反演中，此时待求参数包含电导率和系统飞行姿态角度，
其中，为第1条测线第1个测点第1层的电导率，为第1条测线第1个测点的发射线圈俯仰姿态角度，为第1条测线第1个测点的接收线圈俯仰姿态角度，为第1条测线第1个测点的吊舱同向摆动角度。
空间约束整体反演算法的目标函数Φ(m)由数据误差项Φd(m)，模型先验信息误差项Φr(m)，模型纵向粗糙度Φv(m)，模型横向粗糙度Φh(m)和模型空间粗糙度Φl(m)构成：
Φ(m)＝Φd(m)+λ[1*Φr(m)+50*Φv(mσ)+200*Φh(m)+20*Φl(m)] (7)
λ为正则化因子，在反演迭代过程中通过线性搜索自适应迭代自动确定。
模型先验信息误差项mr为测区航空电磁数据对应的模型先验信息，如步骤a所述，本文采用电导率深度成像结果作为模型先验信息，为先验信息不确定度的方差矩阵。
模型纵向粗糙度由于模型纵向粗糙度仅对单个测点的各层电导率之间进行约束，此时只取m中电导率部分记为mσ，
Qv为模型纵向粗糙度矩阵，对单测点的各层电导率进行约束，本文根据该测点相邻层电导率的二阶差分计算，即
可得模型纵向粗糙度矩阵为，
模型横向粗糙度采用二阶导数计算模型参数的横向粗糙度，对相邻测点间的电导率和飞行姿态角度进行约束，则横向粗糙度矩阵为，
模型横向粗糙度保证了单测线上相邻测点间的连续性，使反演结果沿测线平滑。
模型空间粗糙度同样采用二阶导数计算模型参数的空间粗糙度，对相邻测线间的电导率和飞行姿态角度进行约束，则横向粗糙度矩阵为，
模型空间粗糙度确保了相邻测线间的连续性和反演结果测线间的平滑性。
步骤六所述利用超松弛迭代预处理共轭梯度算法求解迭代方程，得到新反演模型，令目标函数对待求反演模型参数的偏导等于零，求解目标函数的最优解，在第n次迭代过程为，
其中，数据误差项对待求反演模型参数的偏导为：
Gn为雅克比矩阵，为第i个测点对应反演模型的正演结果fforward(mn)对第j个反演模型参数的偏数；模型先验信息误差项对待求反演模型参数的偏导为模型纵向粗糙度对待求反演模型参数的偏导为模型横向粗糙度对待求反演模型参数的偏导为模型空间粗糙度对待求反演模型参数的偏导为整合式(10)为Amn＝b的形式，
利用超松弛迭代预处理共轭梯度算法求解式(13)线性方程组，得到第n次迭代的反演模型参数mn。
步骤七所述的利用新反演模型计算数据误差项，并判断是否达到反演终止条件，利用式(4)计算步骤六得到的反演模型参数mn对应的数据误差项，并带入式(5)，判断是否达到反演终止条件。如满足反演终止条件，进入步骤八，否，回到步骤五，利用步骤六得到的反演模型参数mn建立关于测区测点电导率和姿态角度的拟三维空间约束整体反演目标函数。
步骤八输出反演结果。分别进行拟二维整体反演和拟三维整体反演，反演参数同上，反演结果如图3和图4所示。从图3可以看出，拟二维整体反演结果受噪声干扰严重，反演结果出现多处不连续情况，无法从反演结果中准确分辨异常埋深、边界和电导率。拟三维整体反演增加了各测线间的约束，有效地抑制噪声，从图4可以看到，异常形态清晰，与图2的理论模型一致。

Claims (1)

1.一种固定翼时间域航空电磁数据拟三维空间约束整体反演方法，其特征在于，包括以下步骤：
a、录入固定翼时间域航空电磁数据对应的先验信息；
b、设置待反演参数，包括地下各层电导率和飞行姿态角度初始模型；
c、对各测点数据进行考虑发射线圈俯仰姿态、接收线圈俯仰姿态和吊舱同向摆动的正演计算，在飞机飞行过程中，吊舱的摆动会影响系统的几何参数，发射线圈旋转存在姿态角度时，接收线圈旋转存在姿态角度时，为得到时域电磁响应，首先将频域二次场变换到s域，再通过拉普拉斯逆变换，得到时域感应电动势；
d、利用正演结果计算数据误差项，判断是否达到反演终止条件，是，输出反演结果；否，进行下一步；
e、建立关于测区测点电导率和飞行姿态角度的拟三维空间约束整体反演目标函数，将发射线圈的俯仰姿态角度、接收线圈的俯仰姿态角度和吊舱的同向摆动角度引入拟三维空间约束整体反演中，此时待求参数包含地下各层电导率和系统飞行姿态角度；
f、利用超松弛迭代预处理共轭梯度算法求解迭代方程，得到新反演模型，令目标函数对待求反演模型参数的偏导等于零，利用超松弛迭代预处理共轭梯度算法求解目标函数的最优解；
g、利用新反演模型计算数据误差项，并判断是否达到反演终止条件，否，返回到步骤e，建立关于测区测点电导率和姿态角度的新的拟三维空间约束整体反演目标函数；
h、是，输出反演结果。
CN201810492825.8A 2018-05-22 2018-05-22 固定翼时间域航空电磁数据拟三维空间约束整体反演方法 Pending CN108984818A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810492825.8A CN108984818A (zh) 2018-05-22 2018-05-22 固定翼时间域航空电磁数据拟三维空间约束整体反演方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810492825.8A CN108984818A (zh) 2018-05-22 2018-05-22 固定翼时间域航空电磁数据拟三维空间约束整体反演方法

Publications (1)

Publication Number Publication Date
CN108984818A true CN108984818A (zh) 2018-12-11

Family

ID=64542074

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810492825.8A Pending CN108984818A (zh) 2018-05-22 2018-05-22 固定翼时间域航空电磁数据拟三维空间约束整体反演方法

Country Status (1)

Country Link
CN (1) CN108984818A (zh)

Cited By (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109738715A (zh) * 2019-01-16 2019-05-10 吉林大学 一种磁共振测深频段空间电磁噪声采集装置及方法
CN110244367A (zh) * 2019-06-17 2019-09-17 吉林大学 一种基于地面多基站的ztem系统姿态补偿方法
CN110598367A (zh) * 2019-10-12 2019-12-20 中南大学 一种足迹引导的高效航空电磁法数值模拟方法
CN111126591A (zh) * 2019-10-11 2020-05-08 重庆大学 一种基于空间约束技术的大地电磁深度神经网络反演方法
CN111257954A (zh) * 2020-02-27 2020-06-09 长安大学 一种基于特征反演的车载阵列式探测方法及系统
CN112253090A (zh) * 2020-10-14 2021-01-22 中海油田服务股份有限公司 一种多频电成像的数据参数反演方法和装置
CN112598310A (zh) * 2020-12-29 2021-04-02 北京筹策科技有限公司 一种航班舱位的实时动态分配方法、装置及计算机设备
CN114966870A (zh) * 2022-06-01 2022-08-30 广东省地质物探工程勘察院 一种地面任意回线任意位置多分量联合探测瞬变电磁方法
CN116047614A (zh) * 2022-12-20 2023-05-02 成都理工大学 基于模型空间约束的半航空瞬变电磁数据正则化牛顿反演方法

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20120195165A1 (en) * 2011-01-31 2012-08-02 Chevron U.S.A. Inc. Exploitation of self-consistency and differences between volume images and interpreted spatial/volumetric context
CN102798898A (zh) * 2012-08-20 2012-11-28 中国地质科学院矿产资源研究所 大地电磁场非线性共轭梯度三维反演方法
CN103675927A (zh) * 2013-12-20 2014-03-26 吉林大学 固定翼航空电磁系统接收吊舱摆动角度的校正方法
CN104537714A (zh) * 2015-01-07 2015-04-22 吉林大学 磁共振与瞬变电磁空间约束联合反演方法
CN106199742A (zh) * 2016-06-29 2016-12-07 吉林大学 一种频率域航空电磁法2.5维带地形反演方法
CN107305600A (zh) * 2016-04-21 2017-10-31 新疆维吾尔自治区煤炭科学研究所 最小二乘法电阻率三维近似反演技术

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20120195165A1 (en) * 2011-01-31 2012-08-02 Chevron U.S.A. Inc. Exploitation of self-consistency and differences between volume images and interpreted spatial/volumetric context
CN102798898A (zh) * 2012-08-20 2012-11-28 中国地质科学院矿产资源研究所 大地电磁场非线性共轭梯度三维反演方法
CN103675927A (zh) * 2013-12-20 2014-03-26 吉林大学 固定翼航空电磁系统接收吊舱摆动角度的校正方法
CN104537714A (zh) * 2015-01-07 2015-04-22 吉林大学 磁共振与瞬变电磁空间约束联合反演方法
CN107305600A (zh) * 2016-04-21 2017-10-31 新疆维吾尔自治区煤炭科学研究所 最小二乘法电阻率三维近似反演技术
CN106199742A (zh) * 2016-06-29 2016-12-07 吉林大学 一种频率域航空电磁法2.5维带地形反演方法

Non-Patent Citations (5)

* Cited by examiner, † Cited by third party
Title
LEIF H. COX等: "3D inversion of airborne electromagnetic data using a moving footprint", 《EXPLORATION GEOPHYSICS》 *
QIONGZHANG等: "Airborne electromagnetic data levelling using principal component analysis based on flight line difference", 《JOURNALOFAPPLIEDGEOPHYSICS》 *
刘云鹤等: "三维频率域航空电磁反演研究", 《地球物理学报》 *
王昊等: "固定翼时间域航空电磁系统线圈姿态及吊舱摆动的校正", 《中国地球科学联合学术年会2017》 *
王琦: "固定翼时间域航空电磁数据整体反演", 《中国优秀博硕士学位论文全文数据库（博士） 基础科学辑》 *

Cited By (14)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109738715B (zh) * 2019-01-16 2020-10-02 吉林大学 一种磁共振测深频段空间电磁噪声采集装置及方法
CN109738715A (zh) * 2019-01-16 2019-05-10 吉林大学 一种磁共振测深频段空间电磁噪声采集装置及方法
CN110244367B (zh) * 2019-06-17 2020-05-29 吉林大学 一种基于地面多基站的ztem系统姿态补偿方法
CN110244367A (zh) * 2019-06-17 2019-09-17 吉林大学 一种基于地面多基站的ztem系统姿态补偿方法
CN111126591A (zh) * 2019-10-11 2020-05-08 重庆大学 一种基于空间约束技术的大地电磁深度神经网络反演方法
CN111126591B (zh) * 2019-10-11 2023-04-18 重庆大学 一种基于空间约束技术的大地电磁深度神经网络反演方法
CN110598367A (zh) * 2019-10-12 2019-12-20 中南大学 一种足迹引导的高效航空电磁法数值模拟方法
CN111257954A (zh) * 2020-02-27 2020-06-09 长安大学 一种基于特征反演的车载阵列式探测方法及系统
CN112253090A (zh) * 2020-10-14 2021-01-22 中海油田服务股份有限公司 一种多频电成像的数据参数反演方法和装置
CN112598310A (zh) * 2020-12-29 2021-04-02 北京筹策科技有限公司 一种航班舱位的实时动态分配方法、装置及计算机设备
CN114966870A (zh) * 2022-06-01 2022-08-30 广东省地质物探工程勘察院 一种地面任意回线任意位置多分量联合探测瞬变电磁方法
CN114966870B (zh) * 2022-06-01 2024-03-12 广东省地质物探工程勘察院 一种地面任意回线任意位置多分量联合探测瞬变电磁方法
CN116047614A (zh) * 2022-12-20 2023-05-02 成都理工大学 基于模型空间约束的半航空瞬变电磁数据正则化牛顿反演方法
CN116047614B (zh) * 2022-12-20 2023-10-24 成都理工大学 基于模型空间约束的半航空瞬变电磁数据正则化牛顿反演方法

Similar Documents

Publication Publication Date Title
CN108984818A (zh) 固定翼时间域航空电磁数据拟三维空间约束整体反演方法
CN110058317B (zh) 航空瞬变电磁数据和航空大地电磁数据联合反演方法
CN104537714B (zh) 磁共振与瞬变电磁空间约束联合反演方法
CN108072892B (zh) 一种自动化的地质构造约束层析反演方法
CN108680964A (zh) 一种基于结构约束的归一化重磁电震联合反演方法
CN106772577B (zh) 基于微地震数据和spsa优化算法的震源反演方法
CN106932819B (zh) 基于各向异性马尔科夫随机域的叠前地震参数反演方法
CN110007357B (zh) 一种航空tem和航空mt联合反演方法
CN102393183A (zh) 基于控制网的海量点云快速配准方法
CN106291725B (zh) 一种快速反演地下地质体空间位置的方法
CN104375195A (zh) 时频电磁的多源多分量三维联合反演方法
CN109884710B (zh) 针对激发井深设计的微测井层析成像方法
CN106980736A (zh) 一种各向异性介质的海洋可控源电磁法有限元正演方法
CN114460654B (zh) 基于l1l2混合范数的半航空瞬变电磁数据反演方法及装置
Huang et al. Three-dimensional GPR ray tracing based on wavefront expansion with irregular cells
CN110187382B (zh) 一种回折波和反射波波动方程旅行时反演方法
CN113671570B (zh) 一种地震面波走时和重力异常联合反演方法与系统
WO2012001388A2 (en) Gravity survey data processing
CN104360396B (zh) 一种海上井间tti介质三种初至波走时层析成像方法
CN106483559A (zh) 一种地下速度模型的构建方法
CN108845358B (zh) 断层及构造异常体识别方法及装置
CN109636912A (zh) 应用于三维声呐图像重构的四面体剖分有限元插值方法
CN114114438B (zh) 一种回线源地空瞬变电磁数据的拟三维反演方法
CN110244367B (zh) 一种基于地面多基站的ztem系统姿态补偿方法
CN107356972A (zh) 一种各向异性介质的成像方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
WD01 Invention patent application deemed withdrawn after publication

Application publication date: 20181211

WD01 Invention patent application deemed withdrawn after publication