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

摘要

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

3)计算拟合误差，如果达到设定精度或最大迭代次数，退出计算，否则继续；

4)用拟正演计算雅克比矩阵，得到模型更新步长；

5)更新模型参数，m^k+1＝m^k+Δm。

本发明技术方案的进一步改进在于：所述步骤1)中，定义的目标函数如式1所示，

$\Phi \left(m\right)=\frac{1}{2}\left|\right|d\left(m\right)-d_{obs}=|{|}^2-\frac{\beta }{2}|\left|W(m-m_{ref})\right|{|}^2---1,$

式1中，d为正演模拟得到的数据向量，d_obs为观测数据向量，m为模型参数向量，m_ref为参考模型向量或先验信息模型向量，W为模型光滑度矩阵，β为正则化参数。

本发明技术方案的进一步改进在于：所述步骤1)中，假定初始模型为m_i，代入式1并用泰勒展开，对展开后的式1线性化，并忽略高阶项可得式2，

$\Phi \left(m\right)=\frac{1}{2}\left|\right|(d+\frac{\partial d}{\partial m_i}·\Delta m)-d_{obs}|{|}^2+\frac{\beta }{2}|\left|W(m_i-m_{ref})\right|{|}^2---2,$

式2中，Δm为m-m_i；

对式2求导并令其等于零，得到式3所示的高斯牛顿法模型更新迭代公式，

式3中，H为近似海森矩阵，g为目标函数的梯度，J为雅克比矩阵或灵敏度矩阵，所述灵敏度矩阵的元素表示如4式，

本发明技术方案的进一步改进在于：所述步骤2)的正演计算中，假设时谐因子为e^iωt，将电磁场分解为一次场和二次场，基于二次场的双旋度电场方程表示为式5，

$▿\times ▿\times E_s+{\mathrm{i\omega \mu }}_0{\mathrm{\sigma E}}_s=-{\mathrm{i\omega \mu }}_0{\mathrm{\Delta \sigma E}}_p---5,$

式5中，E_s为二次电场，ω为角频率，μ₀为真空中的磁导率，σ为电导率，E_p为背景场，Δσ为总电导率与背景电导率之差，表示为Δσ＝σ-σ_p，σ_p为背景电导率。

本发明技术方案的进一步改进在于：所述步骤2)中，采用伽辽金加权余量法对式5计算，得式6，

$R_k={\int }_{\Omega }{N}_k·[▿\times ▿\times E_s+{\mathrm{i\omega \mu \sigma E}}_s+{\mathrm{i\omega \mu \Delta \sigma E}}_p]dV=0---6,$

由于式中Ω^e代表一个离散单元，Ne为离散单元数，将式6写成式7所示的离散形式，

本发明技术方案的进一步改进在于：所述式7中K^e，M^e为单元刚度矩阵，表示如式8、式9所示，

$K_{kl}^e={\int }_{{\Omega }^e}(▿\times N_k^e)·(▿\times N_l^e)dV---8,$

$M_{kl}^e={\int }_{{\Omega }^e}{N}_k^e·N_l^{e}dV---9,$

式8、9中为矢量基函数，式8和式9采用27点高斯积分进行计算；

将单元刚度矩阵分配到全局刚度矩阵，得到式10所示的大型线性方程组，

KF＝b 10，

式10中，K为稀疏复对称矩阵，F为场值，b为源项；

采用简单的狄利克雷边界条件，认为在离异常体足够远的边界处二次异常场已经衰减为零，即在边界处满足式11所示，

$n\times E{|}_{\partial \Omega }=0---11,$

求解得到电场后，磁场用法拉第定律求得，即式12所示，

$H={(-{\mathrm{i\omega \mu }}_0)}^{-1}▿\times E---12.$

本发明技术方案的进一步改进在于：所述步骤4)中，采用求解伴随方程的方法来求解雅克比矩阵，其中磁场的计算公式为式13所示，

$H={(-{\mathrm{i\omega \mu }}_0)}^{-1}\underset{i=1}{\overset{12}{\Sigma }}▿\times N_i{E}_i---13,$

磁场对模型参数的导数为式14所示，

$\frac{\partial H}{\partial m}={(-{\mathrm{i\omega \mu }}_0)}^{-1}\underset{i=1}{\overset{12}{\Sigma }}▿\times N_i\frac{\partial E_i}{\partial m}---14,$

磁场对模型参数的导数可以转换成电场对模型参数的导数，由于正演最后得到式15所示的大型复数线性方程组，

Ax＝b 15,

对式(14)两边同时对m求导数，得式16，

$A\frac{\partial x}{\partial m}=-\frac{\partial A}{\partial m}x+\frac{\partial b}{\partial m}---16,$

通过求解式16得到电场对模型参数的导数。

本发明技术方案的进一步改进在于：所述步骤5)中，通过求解线性方程组得到模型参数更新量，即采用伴随方程的方法来求解雅克比矩阵，进而对模型更新方程组求解。

本发明技术方案的进一步改进在于：所述步骤5)中，采用式17更新下一次迭代的模型，

m_i+1＝m_i+αΔm>

式17中，Δm为模型更新向量，α为步长，取值范围为0<α≤1；

采用最速下降公式式18来选取合适的步长，

$\phi \left(m_{i+1}\right)=\phi (m_i+\alpha \Delta m)\le \phi \left(m_i\right)+c_1\alpha ▿\phi \left(m_i\right)\Delta m---18,$

式中c₁为一常量，通常取值为10^-4，开始α取值为1，判断其是否满足式18，若满足则取该α值，更新模型，否则α减少为原来的1/2，重复上述步骤，直到满足为止。

本发明技术方案的进一步改进在于：所述步骤5)中，采用下列方法之一选择模型正则化因子，

①在整个反演过程中β固定为一个值；

②在每次迭代过程中逐渐减小β值。

由于采用上述技术方案，本发明所产生的有益效果在于：

本申请的技术方案考虑带地形反演，消除了地形的影响，实现了起伏地表条件下的频率域航空电磁法2.5维反演，为山区航空电磁数据处理与解释提供了一种行之有效的计算方法。

本发明公开的一种频率域航空电磁法2.5维带地形反演方法。为简化2.5维正反演计算的复杂性，利用三维正演算法来进行二维模型数值模拟，减少了波数域的正反向变换，保证了计算速度与精度。正演利用三维矢量有限元方法进行计算，采用非规则六面体网格进行剖分,确保对起伏地表的精确模拟。

为提高多源电磁法的正反演速度，采用大规模稀疏矩阵并行直接求解器PARDISO来同时求解大型线性方程组及雅克比矩阵。反演方法使用高斯牛顿法，首先采用互换定理来快速求解雅克比矩阵,然后采用共轭梯度法来求解线性方程组，实现了方程的快速求解。

本发明采用三维矢量有限元法进行正演，满足了不同物性边界处电磁场切向连续，可以避免传统节点有限方法出现伪解的不足；同时矢量基函数自动满足散度为零的条件，不需要进行散度校正。有限元方法可以较好的模拟地表起伏，克服有限差分方法在模拟地形时的台阶效应。此外，与传统节点有限元相比，可以减少未知数个数，从而减少线性方程组的维度，一定程度上减少计算时间。针对航空电磁的多源问题，引入直接求解方程组的方法，只需一次矩阵分解，然后只需简单的向前向后迭代。反演采用高斯牛顿反演方案，可以加速收敛，减少迭代次数。通过对比带地形和不带地形反演对比研究，发现地形对航空电磁法的观测结果影响较大，忽略地形往往得不到正确的结果。

所以，本发明实现了频率域航空电磁法的2.5维带地形反演，解决了忽略地形反演造成许多虚假异常的问题，消除了地形因素对航空电磁数据处理的影响，提高了航空电磁法数据处理的准确性，为山区航空电磁数据处理与解释提供了一种行之有效的计算方法。

附图说明

图1单元全局坐标映射到局部坐标图，图中E_i为单元刚度矩阵，单元刚度矩阵按图中规则分配到全局刚度矩阵，从而得到大型线性方程组；

图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得到二次磁场H_x和H_z；

3)计算拟合误差，如果达到设定精度或最大迭代次数，退出计算，否则继续；

4)用拟正演计算雅克比矩阵，得到模型更新步长；

5)更新模型参数，m^k+1＝m^k+Δm。

在步骤1)中，定义的目标函数如式1所示，

$\Phi \left(m\right)=\frac{1}{2}\left|\right|d\left(m\right)-d_{obs}=|{|}^2-\frac{\beta }{2}|\left|W(m-m_{ref})\right|{|}^2---1,$

式1中，d为正演模拟得到的数据向量，d_obs为观测数据向量，m为模型参数向量，m_ref为参考模型向量或先验信息模型向量，W为模型光滑度矩阵，β为正则化参数，β起到均衡数据拟合项和模型约束项的作用。

用式1定义的目标函数由两部分组成，第一部分为数据拟合项，第二部分为正则化项。第一部分确保反演出来的模型能拟合上观测数据，第二部分约束反演结果与已知先验信息之间的差距。采用m-m_ref的二阶差分算子最小的第二部分的光滑模型正则化第二项，即W为二阶差分算子。

假定初始模型为m_i，代入式1并用泰勒展开，对展开后的式1线性化，并忽略高阶项可得式2，

$\Phi \left(m\right)=\frac{1}{2}\left|\right|(d+\frac{\partial d}{\partial m_i}·\Delta m)-d_{obs}|{|}^2+\frac{\beta }{2}|\left|W(m_i-m_{ref})\right|{|}^2---2,$

式2中，Δm为m-m_i；

对式2求导并令其等于零，得到式3所示的高斯牛顿法模型更新迭代公式，

式3中，H为近似海森矩阵，g为目标函数的梯度，J为雅克比矩阵或灵敏度矩阵，所述灵敏度矩阵的元素表示如4式，

在步骤2)的正演计算中，为了消除源的奇异性，假设时谐因子为e^iωt，将电磁场分解为一次场和二次场，一次场即为背景场，二次场即为散射场，基于二次场的双旋度电场方程表示为式5，

$▿\times ▿\times E_s+{\mathrm{i\omega \mu }}_0{\mathrm{\sigma E}}_s=-{\mathrm{i\omega \mu }}_0{\mathrm{\Delta \sigma E}}_p---5,$

式5中，E_s为二次电场，ω为角频率，μ₀为真空中的磁导率，σ为电导率，E_p为背景场，Δσ为总电导率与背景电导率之差，表示为Δσ＝σ-σ_p，σ_p为背景电导率。背景场在均匀全空间或层状介质下用解析解进行计算，二次电场用矢量有限元进行离散计算。

采用伽辽金加权余量法对式5计算，得式6，

$R_k={\int }_{\Omega }{N}_k·[▿\times ▿\times E_s+{\mathrm{i\omega \mu \sigma E}}_s+{\mathrm{i\omega \mu \Delta \sigma E}}_p]dV=0---6,$

由于式中Ω^e代表一个离散单元，Ne为离散单元数，将式6写成式7所示的离散形式，

式7中K^e，M^e为单元刚度矩阵，表示如式8、式9所示，

$K_{kl}^e={\int }_{{\Omega }^e}(▿\times N_k^e)·(▿\times N_l^e)dV---8,$

$M_{kl}^e={\int }_{{\Omega }^e}{N}_k^e·N_l^{e}dV---9,$

式8、9中为矢量基函数，式8和式9采用27点高斯积分进行计算；

将单元刚度矩阵按一定规则分配到全局刚度矩阵(具体规则见附图1)，得到式10所示的大型线性方程组，

KF＝b 10，

式10中，K为稀疏复对称矩阵，F为场值，b为源项；

为了确保解的唯一性，采用简单的狄利克雷边界条件，认为在离异常体足够远的边界处二次异常场已经衰减为零，即在边界处满足式11所示，

$n\times E{|}_{\partial \Omega }=0---11,$

一旦求解得到电场后，磁场用法拉第定律求得，即式12所示，

$H={(-{\mathrm{i\omega \mu }}_0)}^{-1}▿\times E---12.$

在步骤4)中，由式(3)可知，模型迭代更新的关键是求解式中线性方程组，而解线性方程组的关键是求取雅克比矩阵或灵敏度矩阵J。由此可见，雅克比矩阵的计算是反演方案中至关重要的一步，其求解的效率很大程度上决定了反演的效率。雅克比矩阵的计算方法有解析法、差分方法和伴随方程法。由于高维反演是无法用解析方法求解雅克比矩阵的，而差分方法实现简单，但计算效率太低。所以，采用求解伴随方程的方法来求解雅克比矩阵，其中磁场的计算公式为式13所示，

$H={(-{\mathrm{i\omega \mu }}_0)}^{-1}\underset{i=1}{\overset{12}{\Sigma }}▿\times N_i{E}_i---13,$

磁场对模型参数的导数为式14所示，

$\frac{\partial H}{\partial m}={(-{\mathrm{i\omega \mu }}_0)}^{-1}\underset{i=1}{\overset{12}{\Sigma }}▿\times N_i\frac{\partial E_i}{\partial m}---14,$

即磁场对模型参数的导数可以转换成电场对模型参数的导数，由于正演最后得到式15所示的大型复数线性方程组，

Ax＝b 15,

对式(14)两边同时对m求导数，得式16，

$A\frac{\partial x}{\partial m}=-\frac{\partial A}{\partial m}x+\frac{\partial b}{\partial m}---16,$

通过求解式16得到电场对模型参数的导数。

式16与正演线性方程组具有相同的系数矩阵A，但右端项不同。该求解过程类似于正演过程，称其为‘伴随正演’。依据互换定理，只需把单位源放在接收点处进行一次‘伴随正演’，就可以得到接收点处磁场对模型参数的灵敏度。接收点对地下某一模型参数的灵敏度等于把单位源置于该接收点时与该模型参数相关单元场值的简单加权求和，权值可由式右端项得到。这样只需一次伴随正演就可以得到雅克比矩阵。这里一次伴随正演指的是求解完所有测点和频点。

在步骤5)中，由公式(3)可知，要得到模型参数更新量，就必须要求解线性方程组。可以采用直接解法或迭代方法求解该线性方程组。由(3)式可知，求解该方程要计算J^TJ，如果采用直接解法，势必要显示存储雅克比矩阵J。采用预优共轭梯度迭代法来求解线性方程组，只需要计算矩阵和向量的乘积用迭代法解该线性方程组可以采用存储雅克比矩阵或不存储雅克比矩阵两种方案。如果采用存储方案，不需要每次迭代都重新求解雅克比矩阵J，但对于高维反演来说，内存的需求量是很大的，因为雅克比矩阵是稠密矩阵，需要全部存储。如果不存储雅克比矩阵，求解方程组时，每次迭代都要重新求取雅克比矩阵，需额外进行一次正演和伴随正演计算，但可以避免内存不足的问题。

所以，采用预优共轭梯度迭代法来求解线性方程组，只需要计算矩阵和向量的乘积用迭代法解该线性方程组可以采用存储雅克比矩阵或不存储雅克比矩阵两种方案。通过求解线性方程组得到模型参数更新量，即采用伴随方程的方法来求解雅克比矩阵，进而对模型更新方程组求解。

一旦求解式得到了模型更新向量Δm，

采用式17更新下一次迭代的模型，

m_i+1＝m_i+αΔm>

式17中，Δm为模型更新向量，α为步长，取值范围为0<α≤1；

采用最速下降公式式18来选取合适的步长，

$\phi \left(m_{i+1}\right)=\phi (m_i+\alpha \Delta m)\le \phi \left(m_i\right)+c_1\alpha ▿\phi \left(m_i\right)\Delta m---18,$

式中c₁为一常量，通常取值为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和16k Hz三个频率进行计算。计算结果如图(3)所示，Sasaki结果见参考文献(Yutaka S.,Hiroomi N.2003.Topographic effects in frequency-domain helicopter-borne electromagnetics.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)所示。虽然该反演结果也能反映出低阻和高阻异常体的大致形态，但其和真实的模型相差较大。低阻体异常体在横向上有较大拉伸。此外，跟低阻组合模型一样，会出现许多冗余的虚假异常，特别是在山峰正下方会出现高阻异常体，这也进一步证明了忽略地形影响不可避免存在计算误差，从而更加体现了本发明的价值。