CN106199742B - 一种频率域航空电磁法2.5维带地形反演方法 - Google Patents

一种频率域航空电磁法2.5维带地形反演方法 Download PDF

Info

Publication number
CN106199742B
CN106199742B CN201610495418.3A CN201610495418A CN106199742B CN 106199742 B CN106199742 B CN 106199742B CN 201610495418 A CN201610495418 A CN 201610495418A CN 106199742 B CN106199742 B CN 106199742B
Authority
CN
China
Prior art keywords
formula
model
equation
matrix
field
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.)
Active
Application number
CN201610495418.3A
Other languages
English (en)
Other versions
CN106199742A (zh
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
Hebei Electric Power Design and Research Institute
Original Assignee
Jilin University
Hebei Electric Power Design and Research Institute
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, Hebei Electric Power Design and Research Institute filed Critical Jilin University
Priority to CN201610495418.3A priority Critical patent/CN106199742B/zh
Publication of CN106199742A publication Critical patent/CN106199742A/zh
Application granted granted Critical
Publication of CN106199742B publication Critical patent/CN106199742B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • 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

Landscapes

  • Life Sciences & Earth Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Complex Calculations (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了一种频率域航空电磁法2.5维带地形反演方法，包括以下步骤：1)定义目标函数，设置迭代次数为i＝0、拟合精度及最大迭代次数，输入初始模型及反演数据；2)进行正演计算，解正演方程KF＝b得到二次磁场Hx和Hz；3)计算拟合误差，如果达到设定精度或最大迭代次数，退出计算，否则继续；4)用拟正演计算雅克比矩阵，得到模型更新步长；5)更新模型参数，mk+1＝mk+Δm。本发明考虑带地形反演，消除了地形的影响，实现了起伏地表条件下的频率域航空电磁法2.5维反演，为山区航空电磁数据处理与解释提供了一种行之有效的计算方法。

Description

一种频率域航空电磁法2.5维带地形反演方法
技术领域
本发明属于是航空物探领域，具体涉及一种在地形起伏条件下的频率域航空电磁法的2.5维正反演数据处理系统。
背景技术
频率域航空电磁法作为一种重要的地球物理勘探方法，在矿产勘查、地质填图、地下水资源勘查和环境监测等众多领域得到广泛运用。航空电磁法经常在山区作业，这些地区地形起伏较大，对航空电磁响应有严重影响，忽略地形影响会给航空电磁数据解释造成很大误差，只有带地形反演才能消除地形的影响。所以，开展航空电磁法带地形反演是非常必要的。
目前，频率域航空电磁法的解释仍主要以电导率成像技术及层状介质的一维反演为主。但一维反演方法对高维模型的重构能力差，在地表起伏变化的情况下，即使进行了地形校正，用一维反演方法仍不可避免存在计算误差。
因此许多学者在高维反演方面做了一些研究。Wilson et al.(2006)实现了2.5维反演，并将该方法应用于理论和实测数据。Cox et al.(2010)开展了基于footprint技术的三维反演算法，对实测数据进行反演。Liu(2013)开发了基于有限差分的三维反演方案，对非线性共轭梯度(NLCG)及有限内存拟牛顿(LBFGS)方法做了对比研究，得出LBFGS方法更适合于三维频率域航空电磁法反演(Huang.，2016)。Yi and Sasaki(2015)提出了用航空电磁数据和地面直流数据进行联合反演的方案，认为这样可以更好的还原真实模型，提高反演分辨率。但是，以上方法都是假设地表平坦的情况下进行的，没有考虑地形的影响。
发明内容
本发明要解决的技术问题是提供一种在地形起伏条件下的频率域航空电磁法的2.5维正反演数据处理方法，解决忽略地形反演造成许多虚假异常的问题，消除地形因素对航空电磁数据处理的影响，提高航空电磁法数据处理的准确性。
为解决上述技术问题，本发明所采取的技术方案是：
一种频率域航空电磁法2.5维带地形反演方法，包括以下步骤，
1)定义目标函数，设置迭代次数为i＝0、拟合精度及最大迭代次数，输入初始模型及反演数据；
2)进行正演计算，解正演方程KF＝b得到二次磁场Hx和Hz；
3)计算拟合误差，如果达到设定精度或最大迭代次数，退出计算，否则继续；
4)用拟正演计算雅克比矩阵，得到模型更新步长；
5)更新模型参数，mk+1＝mk+Δm。
本发明技术方案的进一步改进在于：所述步骤1)中，定义的目标函数如式1所示，
式1中，d为正演模拟得到的数据向量，dobs为观测数据向量，m为模型参数向量，mref为参考模型向量或先验信息模型向量，W为模型光滑度矩阵，β为正则化参数。
本发明技术方案的进一步改进在于：所述步骤1)中，假定初始模型为mi，代入式1并用泰勒展开，对展开后的式1线性化，并忽略高阶项可得式2，
式2中，Δm为m-mi；
对式2求导并令其等于零，得到式3所示的高斯牛顿法模型更新迭代公式，
式3中，H为近似海森矩阵，g为目标函数的梯度，J为雅克比矩阵或灵敏度矩阵，所述灵敏度矩阵的元素表示如4式，
本发明技术方案的进一步改进在于：所述步骤2)的正演计算中，假设时谐因子为eiωt，将电磁场分解为一次场和二次场，基于二次场的双旋度电场方程表示为式5，
式5中，Es为二次电场，ω为角频率，μ0为真空中的磁导率，σ为电导率，Ep为背景场，Δσ为总电导率与背景电导率之差，表示为Δσ＝σ-σp，σp为背景电导率。
本发明技术方案的进一步改进在于：所述步骤2)中，采用伽辽金加权余量法对式5计算，得式6，
由于式中Ωe代表一个离散单元，Ne为离散单元数，将式6写成式7所示的离散形式，
本发明技术方案的进一步改进在于：所述式7中Ke，Me为单元刚度矩阵，表示如式8、式9所示，
式8、9中为矢量基函数，式8和式9采用27点高斯积分进行计算；
将单元刚度矩阵分配到全局刚度矩阵，得到式10所示的大型线性方程组，
KF＝b 10，
式10中，K为稀疏复对称矩阵，F为场值，b为源项；
采用简单的狄利克雷边界条件，认为在离异常体足够远的边界处二次异常场已经衰减为零，即在边界处满足式11所示，
求解得到电场后，磁场用法拉第定律求得，即式12所示，
本发明技术方案的进一步改进在于：所述步骤4)中，采用求解伴随方程的方法来求解雅克比矩阵，其中磁场的计算公式为式13所示，
磁场对模型参数的导数为式14所示，
磁场对模型参数的导数可以转换成电场对模型参数的导数，由于正演最后得到式15所示的大型复数线性方程组，
Ax＝b 15,
对式(14)两边同时对m求导数，得式16，
通过求解式16得到电场对模型参数的导数。
本发明技术方案的进一步改进在于：所述步骤5)中，通过求解线性方程组得到模型参数更新量，即采用伴随方程的方法来求解雅克比矩阵，进而对模型更新方程组求解。
本发明技术方案的进一步改进在于：所述步骤5)中，采用式17更新下一次迭代的模型，
mi+1＝mi+αΔm l7，
式17中，Δm为模型更新向量，α为步长，取值范围为0<α≤1；
采用最速下降公式式18来选取合适的步长，
式中c1为一常量，通常取值为10-4，开始α取值为1，判断其是否满足式18，若满足则取该α值，更新模型，否则α减少为原来的1/2，重复上述步骤，直到满足为止。
本发明技术方案的进一步改进在于：所述步骤5)中，采用下列方法之一选择模型正则化因子，
①在整个反演过程中β固定为一个值；
②在每次迭代过程中逐渐减小β值。
由于采用上述技术方案，本发明所产生的有益效果在于：
本申请的技术方案考虑带地形反演，消除了地形的影响，实现了起伏地表条件下的频率域航空电磁法2.5维反演，为山区航空电磁数据处理与解释提供了一种行之有效的计算方法。
本发明公开的一种频率域航空电磁法2.5维带地形反演方法。为简化2.5维正反演计算的复杂性，利用三维正演算法来进行二维模型数值模拟，减少了波数域的正反向变换，保证了计算速度与精度。正演利用三维矢量有限元方法进行计算，采用非规则六面体网格进行剖分,确保对起伏地表的精确模拟。
为提高多源电磁法的正反演速度，采用大规模稀疏矩阵并行直接求解器PARDISO来同时求解大型线性方程组及雅克比矩阵。反演方法使用高斯牛顿法，首先采用互换定理来快速求解雅克比矩阵,然后采用共轭梯度法来求解线性方程组，实现了方程的快速求解。
本发明采用三维矢量有限元法进行正演，满足了不同物性边界处电磁场切向连续，可以避免传统节点有限方法出现伪解的不足；同时矢量基函数自动满足散度为零的条件，不需要进行散度校正。有限元方法可以较好的模拟地表起伏，克服有限差分方法在模拟地形时的台阶效应。此外，与传统节点有限元相比，可以减少未知数个数，从而减少线性方程组的维度，一定程度上减少计算时间。针对航空电磁的多源问题，引入直接求解方程组的方法，只需一次矩阵分解，然后只需简单的向前向后迭代。反演采用高斯牛顿反演方案，可以加速收敛，减少迭代次数。通过对比带地形和不带地形反演对比研究，发现地形对航空电磁法的观测结果影响较大，忽略地形往往得不到正确的结果。
所以，本发明实现了频率域航空电磁法的2.5维带地形反演，解决了忽略地形反演造成许多虚假异常的问题，消除了地形因素对航空电磁数据处理的影响，提高了航空电磁法数据处理的准确性，为山区航空电磁数据处理与解释提供了一种行之有效的计算方法。
附图说明
图1单元全局坐标映射到局部坐标图，图中Ei为单元刚度矩阵，单元刚度矩阵按图中规则分配到全局刚度矩阵，从而得到大型线性方程组；
图2是二维梯形山模型图；
图3是本算法正演结果与Sasaki有限差分结果对比图；
图4是水平地表情况下反演多异常体模型图；
图5是水平地表情况下水平共面装置正演结果实、虚分量图；
图6是水平地表情况下垂直共面装置正演结果实、虚分量图；
图7是水平地表情况下水平共面装置反演结果图
图8是水平地表情况下垂直同轴装置反演结果图
图9是山丘低阻组合模型图；
图10是山丘低阻组合模型下水平共面装置正演结果实、虚分量图；
图11是山丘低阻组合模型下带地形山丘低阻组合模型反演结果图；
图12是山丘低阻组合模型下不带地形山丘低阻组合模型反演结果图；
图13是山丘高低阻组合模型图；
图14是山丘高低阻组合模型下水平共面装置正演结果实、虚分量图；
图15是山丘高低阻组合模型下带地形山丘低阻组合模型反演结果图；
图16是山丘高低阻组合模型下不带地形山丘低阻组合模型反演结果图；
图17是本发明的计算流程示意图。
具体实施方式
实施例1
一种频率域航空电磁法2.5维带地形反演方法，包括以下步骤，即反演流程如下：
1)定义目标函数，设置迭代次数为i＝0、拟合精度及最大迭代次数，输入初始模型及反演数据；
2)进行正演计算，解正演方程KF＝b得到二次磁场Hx和Hz；
3)计算拟合误差，如果达到设定精度或最大迭代次数，退出计算，否则继续；
4)用拟正演计算雅克比矩阵，得到模型更新步长；
5)更新模型参数，mk+1＝mk+Δm。
在步骤1)中，定义的目标函数如式1所示，
式1中，d为正演模拟得到的数据向量，dobs为观测数据向量，m为模型参数向量，mref为参考模型向量或先验信息模型向量，W为模型光滑度矩阵，β为正则化参数，β起到均衡数据拟合项和模型约束项的作用。
用式1定义的目标函数由两部分组成，第一部分为数据拟合项，第二部分为正则化项。第一部分确保反演出来的模型能拟合上观测数据，第二部分约束反演结果与已知先验信息之间的差距。采用m-mref的二阶差分算子最小的第二部分的光滑模型正则化第二项，即W为二阶差分算子。
假定初始模型为mi，代入式1并用泰勒展开，对展开后的式1线性化，并忽略高阶项可得式2，
式2中，Δm为m-mi；
对式2求导并令其等于零，得到式3所示的高斯牛顿法模型更新迭代公式，
式3中，H为近似海森矩阵，g为目标函数的梯度，J为雅克比矩阵或灵敏度矩阵，所述灵敏度矩阵的元素表示如4式，
在步骤2)的正演计算中，为了消除源的奇异性，假设时谐因子为eiωt，将电磁场分解为一次场和二次场，一次场即为背景场，二次场即为散射场，基于二次场的双旋度电场方程表示为式5，
式5中，Es为二次电场，ω为角频率，μ0为真空中的磁导率，σ为电导率，Ep为背景场，Δσ为总电导率与背景电导率之差，表示为Δσ＝σ-σp，σp为背景电导率。背景场在均匀全空间或层状介质下用解析解进行计算，二次电场用矢量有限元进行离散计算。
采用伽辽金加权余量法对式5计算，得式6，
由于式中Ωe代表一个离散单元，Ne为离散单元数，将式6写成式7所示的离散形式，
式7中Ke，Me为单元刚度矩阵，表示如式8、式9所示，
式8、9中为矢量基函数，式8和式9采用27点高斯积分进行计算；
将单元刚度矩阵按一定规则分配到全局刚度矩阵(具体规则见附图1)，得到式10所示的大型线性方程组，
KF＝b 10，
式10中，K为稀疏复对称矩阵，F为场值，b为源项；
为了确保解的唯一性，采用简单的狄利克雷边界条件，认为在离异常体足够远的边界处二次异常场已经衰减为零，即在边界处满足式11所示，
一旦求解得到电场后，磁场用法拉第定律求得，即式12所示，
在步骤4)中，由式(3)可知，模型迭代更新的关键是求解式中线性方程组，而解线性方程组的关键是求取雅克比矩阵或灵敏度矩阵J。由此可见，雅克比矩阵的计算是反演方案中至关重要的一步，其求解的效率很大程度上决定了反演的效率。雅克比矩阵的计算方法有解析法、差分方法和伴随方程法。由于高维反演是无法用解析方法求解雅克比矩阵的，而差分方法实现简单，但计算效率太低。所以，采用求解伴随方程的方法来求解雅克比矩阵，其中磁场的计算公式为式13所示，
磁场对模型参数的导数为式14所示，
即磁场对模型参数的导数可以转换成电场对模型参数的导数，由于正演最后得到式15所示的大型复数线性方程组，
Ax＝b 15,
对式(14)两边同时对m求导数，得式16，
通过求解式16得到电场对模型参数的导数。
式16与正演线性方程组具有相同的系数矩阵A，但右端项不同。该求解过程类似于正演过程，称其为‘伴随正演’。依据互换定理，只需把单位源放在接收点处进行一次‘伴随正演’，就可以得到接收点处磁场对模型参数的灵敏度。接收点对地下某一模型参数的灵敏度等于把单位源置于该接收点时与该模型参数相关单元场值的简单加权求和，权值可由式右端项得到。这样只需一次伴随正演就可以得到雅克比矩阵。这里一次伴随正演指的是求解完所有测点和频点。
在步骤5)中，由公式(3)可知，要得到模型参数更新量，就必须要求解线性方程组。可以采用直接解法或迭代方法求解该线性方程组。由(3)式可知，求解该方程要计算JTJ，如果采用直接解法，势必要显示存储雅克比矩阵J。采用预优共轭梯度迭代法来求解线性方程组，只需要计算矩阵和向量的乘积用迭代法解该线性方程组可以采用存储雅克比矩阵或不存储雅克比矩阵两种方案。如果采用存储方案，不需要每次迭代都重新求解雅克比矩阵J，但对于高维反演来说，内存的需求量是很大的，因为雅克比矩阵是稠密矩阵，需要全部存储。如果不存储雅克比矩阵，求解方程组时，每次迭代都要重新求取雅克比矩阵，需额外进行一次正演和伴随正演计算，但可以避免内存不足的问题。
所以，采用预优共轭梯度迭代法来求解线性方程组，只需要计算矩阵和向量的乘积用迭代法解该线性方程组可以采用存储雅克比矩阵或不存储雅克比矩阵两种方案。通过求解线性方程组得到模型参数更新量，即采用伴随方程的方法来求解雅克比矩阵，进而对模型更新方程组求解。
一旦求解式得到了模型更新向量Δm，
采用式17更新下一次迭代的模型，
mi+1＝mi+αΔm 17，
式17中，Δm为模型更新向量，α为步长，取值范围为0<α≤1；
采用最速下降公式式18来选取合适的步长，
式中c1为一常量，通常取值为10-4，开始α取值为1，判断其是否满足式18，若满足则取该α值，更新模型，否则α减少为原来的1/2，重复上述步骤，直到满足为止。
正则化因子β在反演中起到均衡数据拟合项和模型约束项的权重作用，目前，正则化因子β的选择，主要有3种方案，包括：
①在整个反演过程中β固定为一个值；
②在每次迭代过程中逐渐减小β值。
③在每次迭代中用L曲线法选择最合适的β值。
采用方案①或方案②选择模型正则化因子。
实施例2
步骤2)中进行的正演计算，正演精度验证的过程如下：
为了验证本算法的正确性，采用和Sasaki一样的二维梯形山模型(Yutaka etal.，2003)，如图(2)所示。梯形山的顶部和底部分别为20m和220m，山顶到山底的距离为50m，坡度为26.5°，介质的电阻率为100ohm-m。将模拟区域在x、y和z方向剖分为41×41×41网格，中间区域的网格大小为10m×10m×10m，边界处的最大网格大小为1280m×1280m×1280m。用水平共面装置(HCP)进行计算，收发线圈距地表30m，收发距为10m，用1k Hz、4k Hz和16kHz三个频率进行计算。计算结果如图(3)所示，Sasaki结果见参考文献(Yutaka S.,HiroomiN.2003.Topographic effects in frequency-domain helicopter-borneelectromagnetics.Exploration Geophysics,34(2):24-28.)，计算结果和Sasaki的有限差分的结果高度吻合，表明本算法具有较高的计算精度。
反演算法验证的过程如下：为了验证本算法的有效性，设计一个均匀半空间模型来进行验证试算，将4个低阻异常体埋藏在均匀半空间中如图(4)所示。低阻异常体的电阻率为10ohm-m，背景电阻率为1000ohm-m。采用水平共面(HCP)和垂直同轴(VCX)装置进行正演，共采用四个频率，分别为1000Hz、2700Hz、7400Hz和20000Hz。线圈收发距离8m，离地面高度为20m，点距为10m，共51个测点。将模拟区域剖分为69×34×29网格，中间区网格大小为10m×10m×10m，扩边区最大网格为1280m×1280m×1280m。水平共面(HCP)和垂直同轴(VCX)装置正演结果分别如图(5)和图(6)所示，其值用一次场进行归一化，单位为ppm。
反演模型剖分为52×1×10块，共计520块未知电阻率参数。初始模型为300ohm-m的均匀半空间模型。经10次迭代，HCP收敛为原始的1.8％，VCX收敛为原始的0.7％。水平共面(HCP)装置反演结果如图(7)所示。垂直同轴(VCX)装置反演结果如图(8)所示。从图(7)和图(8)可以看出，无论是HCP还是VCX都较好的重构了真实模型，结果较好。
实施例3
设计一个山丘低阻组合模型，如图(9)所示。有两个低阻异常体埋藏在均匀半空间，异常体的电阻率为10ohm-m，背景电阻率为300ohm-m。山顶宽20m，
距山脚的距离为50m，坡度为26.6°。我们通过改变网格的z轴坐标来模拟起伏地表。采用水平共面(HCP)装置进行正演，共采用四个频率，分别为1000Hz、2700Hz、7400Hz和20000Hz。线圈收发距离8m，离地面高度为20m，点距为10m，共51个测点。将模拟区域剖分为69×34×29网格，中间区网格大小为10m×10m×10m，扩边区最大网格为1280m×1280m×1280m。水平共面(HCP)装置正演结果如图(10)所示，其值用一次场进行归一化，单位为ppm。
反演模型剖分为52×1×10块，共计520块未知电阻率参数。初始模型为300ohm-m的均匀半空间模型。经10次迭代收敛到0.5％。水平共面(HCP)装置反演结果如图(11)所示。从图(11)可以看出，反演结果很好地还原了真实的电阻率模型，反演效果较好。
为了说明带地形反演的有效性，我们给出了不带地形的反演结果，如图(12)所示。从图(12)可以看出，忽略地形的反演结果并不理想。虽然不带地形反演也能反映出低阻异常体，但其位置和真实模型相差较远，特别是在横向上的拉伸。此外，在反演剖面中会出现许多虚假异常，特别是在山峰正下方会出现一个高阻异常体，这是地形影响造成的。
实施例4
设计一个山丘高低阻组合模型，如图(13)所示。低阻异常体电阻率值为10ohm-m，高阻异常体电阻率值为3000ohm-m。其它参数和低阻组合模型完全一样。正演结果如图(14)所示。
反演初始电阻率为300ohm-m均匀半空间模型。经过10次迭代，拟合差从降至第一次迭代的2.2％。带地形反演结果如图(15)所示。从图(15)可以看出，反演结果较好的还原了真实地电模型，没有多余的异常。忽略地形的反演结果如图(16)所示。虽然该反演结果也能反映出低阻和高阻异常体的大致形态，但其和真实的模型相差较大。低阻体异常体在横向上有较大拉伸。此外，跟低阻组合模型一样，会出现许多冗余的虚假异常，特别是在山峰正下方会出现高阻异常体，这也进一步证明了忽略地形影响不可避免存在计算误差，从而更加体现了本发明的价值。

Claims (1)

1.一种频率域航空电磁法2.5维带地形反演方法，其特征在于：包括以下步骤，
1)定义目标函数，设置迭代次数为i＝0、拟合精度及最大迭代次数，输入初始模型及反演数据；
2)进行正演计算，解正演方程KF＝b得到二次磁场Hx和Hz；
3)计算拟合误差，如果达到设定精度或最大迭代次数，退出计算，否则继续；
4)用拟正演计算雅克比矩阵，得到模型更新步长；
5)更新模型参数，mk+1＝mk+Δm；
所述步骤1)中，定义的目标函数如式(1)所示，
<mrow> <mi>&amp;Phi;</mi> <mrow> <mo>(</mo> <mi>m</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> <mo>|</mo> <mo>|</mo> <mi>d</mi> <mrow> <mo>(</mo> <mi>m</mi> <mo>)</mo> </mrow> <mo>-</mo> <msub> <mi>d</mi> <mrow> <mi>o</mi> <mi>b</mi> <mi>s</mi> </mrow> </msub> <mo>|</mo> <msup> <mo>|</mo> <mn>2</mn> </msup> <mo>+</mo> <mfrac> <mi>&amp;beta;</mi> <mn>2</mn> </mfrac> <mo>|</mo> <mo>|</mo> <mi>W</mi> <mrow> <mo>(</mo> <mi>m</mi> <mo>-</mo> <msub> <mi>m</mi> <mrow> <mi>r</mi> <mi>e</mi> <mi>f</mi> </mrow> </msub> <mo>)</mo> </mrow> <mo>|</mo> <msup> <mo>|</mo> <mn>2</mn> </msup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>,</mo> </mrow>
式(1)中，d为正演模拟得到的数据向量，dobs为观测数据向量，m为模型参数向量，mref为参考模型向量或先验信息模型向量，W为模型光滑度矩阵，β为正则化参数；
所述步骤1)中，假定初始模型为mi，代入式(1)并用泰勒展开，对展开后的式(1)线性化，并忽略高阶项可得式(2)，
<mrow> <mi>&amp;Phi;</mi> <mrow> <mo>(</mo> <mi>m</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> <mo>|</mo> <mo>|</mo> <mrow> <mo>(</mo> <mi>d</mi> <mo>+</mo> <mfrac> <mrow> <mo>&amp;part;</mo> <mi>d</mi> </mrow> <mrow> <mo>&amp;part;</mo> <msub> <mi>m</mi> <mi>i</mi> </msub> </mrow> </mfrac> <mo>&amp;CenterDot;</mo> <mi>&amp;Delta;</mi> <mi>m</mi> <mo>)</mo> </mrow> <mo>-</mo> <msub> <mi>d</mi> <mrow> <mi>o</mi> <mi>b</mi> <mi>s</mi> </mrow> </msub> <mo>|</mo> <msup> <mo>|</mo> <mn>2</mn> </msup> <mo>+</mo> <mfrac> <mi>&amp;beta;</mi> <mn>2</mn> </mfrac> <mo>|</mo> <mo>|</mo> <mi>W</mi> <mrow> <mo>(</mo> <msub> <mi>m</mi> <mi>i</mi> </msub> <mo>-</mo> <msub> <mi>m</mi> <mrow> <mi>r</mi> <mi>e</mi> <mi>f</mi> </mrow> </msub> <mo>)</mo> </mrow> <mo>|</mo> <msup> <mo>|</mo> <mn>2</mn> </msup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>2</mn> <mo>)</mo> </mrow> <mo>,</mo> </mrow>
式(2)中，Δm为m-mi；
对式(2)求导并令其等于零，得到式(3)所示的高斯牛顿法模型更新迭代公式，
式(3)中，H为近似海森矩阵，g为目标函数的梯度，J为雅克比矩阵或灵敏度矩阵，所述灵敏度矩阵的元素表示如(4)式，
所述步骤2)的正演计算中，假设时谐因子为eiωt，将电磁场分解为一次场和二次场，基于二次场的双旋度电场方程表示为式(5)，
▽×▽×Es+iωμ0σEs＝-iωμ0ΔσEp (5)，
式(5)中，Es为二次电场，ω为角频率，μ0为真空中的磁导率，σ为电导率，Ep为背景场，Δσ为总电导率与背景电导率之差，表示为Δσ＝σ-σp，σp为背景电导率；
所述步骤2)中，采用伽辽金加权余量法对式(5)计算，得式(6)，
Rk＝∫ΩNk·[▽×▽×Es+iωμσEs+iωμΔσEp]dV＝0 (6)，由于式中Ωe代表一个离散单元，Nk为离散单元数，将式(6)写成式(7)所示的离散形式，
<mrow> <msub> <mi>R</mi> <mi>k</mi> </msub> <mo>=</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>e</mi> <mo>=</mo> <mn>1</mn> </mrow> <msub> <mi>N</mi> <mi>e</mi> </msub> </munderover> <mo>&amp;lsqb;</mo> <msup> <mi>K</mi> <mi>e</mi> </msup> <msubsup> <mi>E</mi> <mi>s</mi> <mi>e</mi> </msubsup> <mo>-</mo> <msup> <mi>i&amp;omega;&amp;mu;&amp;sigma;M</mi> <mi>e</mi> </msup> <msubsup> <mi>E</mi> <mi>s</mi> <mi>e</mi> </msubsup> <mo>-</mo> <msup> <mi>i&amp;omega;&amp;mu;&amp;Delta;&amp;sigma;M</mi> <mi>e</mi> </msup> <msubsup> <mi>E</mi> <mi>p</mi> <mi>e</mi> </msubsup> <mo>&amp;rsqb;</mo> <mo>=</mo> <mn>0</mn> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>7</mn> <mo>)</mo> </mrow> <mo>;</mo> </mrow>
所述式(7)中Ke，Me为单元刚度矩阵，表示如式(8)、式(9)所示，
<mrow> <msubsup> <mi>K</mi> <mrow> <mi>k</mi> <mi>l</mi> </mrow> <mi>e</mi> </msubsup> <mo>=</mo> <msub> <mo>&amp;Integral;</mo> <msup> <mi>&amp;Omega;</mi> <mi>e</mi> </msup> </msub> <mrow> <mo>(</mo> <mo>&amp;dtri;</mo> <mo>&amp;times;</mo> <msubsup> <mi>N</mi> <mi>k</mi> <mi>e</mi> </msubsup> <mo>)</mo> </mrow> <mo>&amp;CenterDot;</mo> <mrow> <mo>(</mo> <mo>&amp;dtri;</mo> <mo>&amp;times;</mo> <msubsup> <mi>N</mi> <mi>l</mi> <mi>e</mi> </msubsup> <mo>)</mo> </mrow> <mi>d</mi> <mi>V</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mo>(</mo> <mn>8</mn> <mo>)</mo> <mo>,</mo> </mrow>
<mrow> <msubsup> <mi>M</mi> <mrow> <mi>k</mi> <mi>l</mi> </mrow> <mi>e</mi> </msubsup> <mo>=</mo> <msub> <mo>&amp;Integral;</mo> <msup> <mi>&amp;Omega;</mi> <mi>e</mi> </msup> </msub> <msubsup> <mi>N</mi> <mi>k</mi> <mi>e</mi> </msubsup> <mo>&amp;CenterDot;</mo> <msubsup> <mi>N</mi> <mi>l</mi> <mi>e</mi> </msubsup> <mi>d</mi> <mi>V</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>9</mn> <mo>)</mo> </mrow> <mo>,</mo> </mrow>
式(8)、(9)中为矢量基函数，式(8)和式(9)采用27点高斯积分进行计算；
将单元刚度矩阵分配到全局刚度矩阵，得到式(10)所示的大型线性方程组，
KF＝b (10)，
式(10)中，K为稀疏复对称矩阵，F为场值，b为源项；
采用简单的狄利克雷边界条件，认为在离异常体足够远的边界处二次异常场已经衰减为零，即在边界处满足式(11)所示，
<mrow> <mi>n</mi> <mo>&amp;times;</mo> <mi>E</mi> <msub> <mo>|</mo> <mrow> <mo>&amp;part;</mo> <mi>&amp;Omega;</mi> </mrow> </msub> <mo>=</mo> <mn>0</mn> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>11</mn> <mo>)</mo> </mrow> <mo>,</mo> </mrow>
求解得到电场后，磁场用法拉第定律求得，即式(12)所示，
H＝(-iωμ0)-1▽×E (12)；
所述步骤4)中，采用求解伴随方程的方法来求解雅克比矩阵，其中磁场的计算公式为式(13)所示，
<mrow> <mi>H</mi> <mo>=</mo> <msup> <mrow> <mo>(</mo> <mo>-</mo> <msub> <mi>i&amp;omega;&amp;mu;</mi> <mn>0</mn> </msub> <mo>)</mo> </mrow> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msup> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mn>12</mn> </munderover> <mo>&amp;dtri;</mo> <mo>&amp;times;</mo> <msubsup> <mi>N</mi> <mi>i</mi> <mi>e</mi> </msubsup> <msub> <mi>E</mi> <mi>i</mi> </msub> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>13</mn> <mo>)</mo> </mrow> <mo>,</mo> </mrow>
磁场对模型参数的导数为式(14)所示，
<mrow> <mfrac> <mrow> <mo>&amp;part;</mo> <mi>H</mi> </mrow> <mrow> <mo>&amp;part;</mo> <mi>m</mi> </mrow> </mfrac> <mo>=</mo> <msup> <mrow> <mo>(</mo> <mo>-</mo> <msub> <mi>i&amp;omega;&amp;mu;</mi> <mn>0</mn> </msub> <mo>)</mo> </mrow> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msup> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mn>12</mn> </munderover> <mo>&amp;dtri;</mo> <mo>&amp;times;</mo> <msub> <mi>N</mi> <mi>i</mi> </msub> <mfrac> <mrow> <mo>&amp;part;</mo> <msub> <mi>E</mi> <mi>i</mi> </msub> </mrow> <mrow> <mo>&amp;part;</mo> <mi>m</mi> </mrow> </mfrac> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>14</mn> <mo>)</mo> </mrow> <mo>,</mo> </mrow>
磁场对模型参数的导数可以转换成电场对模型参数的导数，由于正演最后得到式(15)所示的大型复数线性方程组，
KF＝b (15)，
对式(14)两边同时对m求导数，得式(16)，
<mrow> <mi>K</mi> <mfrac> <mrow> <mo>&amp;part;</mo> <mi>F</mi> </mrow> <mrow> <mo>&amp;part;</mo> <mi>m</mi> </mrow> </mfrac> <mo>=</mo> <mo>-</mo> <mfrac> <mrow> <mo>&amp;part;</mo> <mi>K</mi> </mrow> <mrow> <mo>&amp;part;</mo> <mi>m</mi> </mrow> </mfrac> <mi>F</mi> <mo>+</mo> <mfrac> <mrow> <mo>&amp;part;</mo> <mi>b</mi> </mrow> <mrow> <mo>&amp;part;</mo> <mi>m</mi> </mrow> </mfrac> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>16</mn> <mo>)</mo> </mrow> <mo>,</mo> </mrow>
通过求解式(16)得到电场对模型参数的导数；
所述步骤5)中，通过求解线性方程组得到模型参数更新量，即采用伴随方程的方法来求解雅克比矩阵，进而对模型更新方程组求解；
所述步骤5)中，采用式(17)更新下一次迭代的模型，
mi+1＝mi+αΔm (17)，
式(17)中，Δm为模型更新向量，α为步长，取值范围为0<α≤1；
采用最速下降公式式(18)来选取合适的步长，
φ(mi+1)＝φ(mi+αΔm)≤φ(mi)+c1α▽φ(mi)Δm (18)，
式中c1为一常量，通常取值为10-4，开始α取值为1，判断其是否满足式(18)，若满足则取该α值，更新模型，否则α减少为原来的1/2，重复上述步骤，直到满足为止；
所述步骤5)中，采用下列方法之一选择模型正则化因子，
①在整个反演过程中β固定为一个值；
②在每次迭代过程中逐渐减小β值。
CN201610495418.3A 2016-06-29 2016-06-29 一种频率域航空电磁法2.5维带地形反演方法 Active CN106199742B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201610495418.3A CN106199742B (zh) 2016-06-29 2016-06-29 一种频率域航空电磁法2.5维带地形反演方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201610495418.3A CN106199742B (zh) 2016-06-29 2016-06-29 一种频率域航空电磁法2.5维带地形反演方法

Publications (2)

Publication Number Publication Date
CN106199742A CN106199742A (zh) 2016-12-07
CN106199742B true CN106199742B (zh) 2018-02-02

Family

ID=57462438

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201610495418.3A Active CN106199742B (zh) 2016-06-29 2016-06-29 一种频率域航空电磁法2.5维带地形反演方法

Country Status (1)

Country Link
CN (1) CN106199742B (zh)

Families Citing this family (14)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107121706A (zh) * 2017-05-08 2017-09-01 厦门大学 基于波恩迭代法的航空瞬变电磁电导率三维反演方法
CN107797148B (zh) * 2017-10-18 2018-08-21 中国国土资源航空物探遥感中心 一种基于三维地质建模的航磁异常场分离方法及系统
CN108873084B (zh) * 2018-05-10 2019-10-08 中南大学 一种基于单位分解积分的直流电阻率无单元正演方法
CN108984818A (zh) * 2018-05-22 2018-12-11 吉林大学 固定翼时间域航空电磁数据拟三维空间约束整体反演方法
CN108761397B (zh) * 2018-05-30 2022-05-27 中南大学 基于电磁散射模拟的极化sar模型分解评价方法
CN109238021B (zh) * 2018-09-30 2023-05-26 中国地质大学（武汉） 基于多频坐标轴旋转的地雷探测背景干扰抑制方法及系统
CN110598367A (zh) * 2019-10-12 2019-12-20 中南大学 一种足迹引导的高效航空电磁法数值模拟方法
CN115407423B (zh) * 2021-05-27 2024-02-06 中国自然资源航空物探遥感中心 重、磁测数据三维反演方法及装置
CN114488314B (zh) * 2021-12-10 2022-12-23 成都理工大学 一种基于陆地与水下直流电联合测量的地质反演方法
CN114779355B (zh) * 2022-02-24 2024-04-16 中国地质大学（武汉） 基于发射电流全波形的地面瞬变电磁法反演方法及装置
CN114297905B (zh) * 2022-03-10 2022-06-03 中南大学 一种二维大地电磁场的快速数值模拟方法
CN115113286B (zh) * 2022-07-06 2023-09-15 长江大学 基于多分量频率域航空电磁数据融合三维反演方法
CN115542408B (zh) * 2022-12-05 2023-03-28 成都理工大学 一种海洋瞬变电磁数据预处理与快速定区成像方法
CN116908928B (zh) * 2023-05-15 2024-03-26 重庆大学 一种基于地层自适应加密的大地电磁反演方法

Family Cites Families (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN100487494C (zh) * 2006-01-26 2009-05-13 中国石油天然气集团公司 三维电阻率级联自动成像方法
CN101520517B (zh) * 2008-02-25 2011-06-22 中国石油集团东方地球物理勘探有限责任公司 一种能准确评价碎屑岩盆地含油气目标的方法
WO2014077721A1 (en) * 2012-11-15 2014-05-22 Baker Hughes Incorporated Apparatus and method for downhole transient resistivity measurement and inversion
CN103135140B (zh) * 2013-01-31 2015-11-18 江苏大学 一种无边缘效应的中心回线tem全期真电阻率计算方法
CN104597506B (zh) * 2015-01-26 2017-03-29 吉林大学 频率域地空电磁勘探方法
CN105549099A (zh) * 2015-12-11 2016-05-04 中国石油大学(华东) 基于全空间正则化下延数据的视磁化强度三维反演方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
2.5维起伏地表条件下时间域航空电磁正演模拟;殷长春 等;《地球物理学报》;20150415;第58卷(第4期);第1411-1424页 *
王卫平 等.频率域航空电磁法地形影响和校正方法.《吉林大学学报（地球科学版）》.2015,第45卷(第3期), *
频率域可控源电磁法2.5D正反演;戴世坤 等;《中国有色金属学报》;20130915;第23卷(第9期);第2513-2523页 *

Also Published As

Publication number Publication date
CN106199742A (zh) 2016-12-07

Similar Documents

Publication Publication Date Title
CN106199742B (zh) 一种频率域航空电磁法2.5维带地形反演方法
CN106980736B (zh) 一种各向异性介质的海洋可控源电磁法有限元正演方法
Kadu et al. Salt reconstruction in full-waveform inversion with a parametric level-set method
CN105549106B (zh) 一种重力多界面反演方法
CN110058315B (zh) 一种三维各向异性射频大地电磁自适应有限元正演方法
CN107657137B (zh) 一种有理函数逼近的分数阶电磁反常扩散三维模拟方法
Martin et al. Seismic modeling with radial-basis-function-generated finite differences
CN112949134B (zh) 基于非结构有限元方法的地-井瞬变电磁反演方法
Yin et al. 3D time-domain airborne EM forward modeling with topography
CN105717547A (zh) 一种各向异性介质大地电磁无网格数值模拟方法
NO343091B1 (no) Elektromagnetisk avbildning av undergrunnen ved fire-dimensjonal parallell databehandling
Zhang et al. 3D inversion of time-domain electromagnetic data using finite elements and a triple mesh formulation
Udosen et al. A finite integration forward solver and a domain search reconstruction solver for electrical resistivity tomography (ERT)
Wang et al. 2D joint inversion of CSAMT and magnetic data based on cross-gradient theory
Wang et al. Divergence correction schemes in finite difference method for 3D tensor CSAMT in axial anisotropic media
Manchuk et al. Implementation aspects of sequential Gaussian simulation on irregular points
CN110826283A (zh) 一种预处理器及基于此预处理器的三维有限差分电磁正演计算方法
Feng et al. Contrast between 2D inversion and 3D inversion based on 2D high-density resistivity data
CN114114438B (zh) 一种回线源地空瞬变电磁数据的拟三维反演方法
Yin et al. A fast 3D gravity forward algorithm based on circular convolution
Cheng et al. 3D Step-by-step inversion strategy for audio magnetotellurics data based on unstructured mesh
Cox et al. 3D airborne electromagnetic inversion using a hybrid edge-based FE-IE method with moving sensitivity domain
Wang et al. [Retracted] Three‐Dimensional Inversion of Borehole‐Surface Resistivity Method Based on the Unstructured Finite Element
Zhou et al. A topography-dependent eikonal solver for accurate and efficient computation of traveltimes and their derivatives in 3D heterogeneous media
Kai et al. Optimization method of first-arrival waveform inversion based on the L-BFGS algorithm

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant