CN105589108B - 基于不同约束条件的瞬变电磁快速三维反演方法 - Google Patents
基于不同约束条件的瞬变电磁快速三维反演方法 Download PDFInfo
- Publication number
- CN105589108B CN105589108B CN201510926452.7A CN201510926452A CN105589108B CN 105589108 B CN105589108 B CN 105589108B CN 201510926452 A CN201510926452 A CN 201510926452A CN 105589108 B CN105589108 B CN 105589108B
- Authority
- CN
- China
- Prior art keywords
- mrow
- msub
- infinitesimal
- time constant
- moment
- 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
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)
- Geophysics And Detection Of Objects (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明提供了一种基于不同约束条件的瞬变电磁快速三维反演方法，采用了基于瞬变电磁矩变换的数据处理及三维正演方法，解决了当前三维反演面临的数据量大、正演复杂的难题；构建了反演中时间常数向量的约束条件，采用最优化算法对反演问题进行迭代，使得反演过程向贴近于实际地下结构的方向优化。本发明提出的方法可在普通计算机上实现地下异常目标的快速三维反演，在实时瞬变电磁数据解释工作中具有重要应用前景。
Description
技术领域
本发明涉及地球物理勘探技术地下地质体探测技术领域，尤其涉及一种基于不同约束条件的瞬变电磁快速三维反演方法。
背景技术
瞬变电磁法利用接地电极或不接地回线向地下发射双极性脉冲电流，地下介质在其激发下的感应涡流产生随时间变化的二次场，在一次场的间歇期，使用接收线圈测量磁场信号，通过对二次场信号的提取和分析，从而达到探测地下地质体的目的。目前作为主要的非地震方法之一，广泛应用于油气、矿产等地下资源探测领域。瞬变电磁数据的反演解释工作是瞬变电磁法勘探中的重要环节。由于高维正演算法的复杂性，多维反演问题尚未妥善解决，实际应用中，对瞬变电磁数据的解释工作主要集中在一维反演。
但是高维反演能够提供更为精细的地电结构信息，目前，随着计算技术的发展，国外已经开展三维反演工作，主要利用积分方程法和有限元法等方法实现严格三维正反演。瞬变电磁的严格三维反演方法受限于复杂的三维正演算法，数据量巨大，占用资源过多，在普通计算机上几乎无法运行，同时严格三维反演的运算速度缓慢，需要几个小时甚至几天才能完成三维反演，由于这些限制，严格三维反演还不能真正投入实际应用，也无法实时进行瞬变电磁数据的解释。因此，本领域技术人员需要迫切解决的一个技术问题就是：如何实现一种可投入应用的三维反演方法，能够高效及准确地反演出地下异常目标体。
发明内容
(一)要解决的技术问题
为了解决现有的三维反演方法中存在的数据量大、运算速度慢、难以投入应用的技术问题，本发明提供了一种基于不同约束条件的瞬变电磁快速三维反演方法。
(二)技术方案
根据本发明的一个方面，提供了一种基于不同约束条件的瞬变电磁快速三维反演方法，包括：步骤A：在地面上铺设发射装置和接收点，设定地下异常区域，将异常区域划分为立方体形状的微元，根据发射装置、接收点以及微元的几何参数计算几何耦合因子矩阵；步骤B：发射装置发射电流信号，电流关断后，各个接收点采集磁场数据，然后采用均匀大地模型，将采集的磁场数据转换为视电导率深度图；步骤C：根据瞬变电磁一阶矩变换，基于视电导率深度图和测量磁场数据获取接收点处的异常区域的参考一阶矩；步骤D：构建反演中时间常数向量的约束条件；步骤E：基于异常区域的参考一阶矩及其误差、微元对应的时间常数向量以及几何耦合因子矩阵构建反演运算的目标函数；步骤F：基于微元对应的时间常数向量的约束条件，采用最优化算法对反演运算的目标函数进行迭代；以及步骤G：判断迭代后的目标函数是否收敛，若收敛，则保存最优化时间常数向量作为反演最终结果，即τ＝(τ1，τ2，…，τK)T，否则，返回步骤F执行，其中，τk是异常区域内第k个微元的电导率对应的时间常数，根据时间常数向量的反演最终结果得到异常体的位置和体积。
(三)有益效果
从上述技术方案可以看出，本发明的基于不同约束条件的瞬变电磁快速三维反演方法具有以下有益效果：
(1)采用了基于瞬变电磁矩变换的数据处理方法，使得三维反演过程处理的数据量大幅降低，加快了运算速度，并且使三维反演能够在普通计算机上运行；
(2)采用了简化的三维正演算法，从而使反演过程中正演部分的运算加快；
(3)采用了基于约束条件的最优化算法，使反演过程向贴近于实际地下结构的方向优化，提高了探测精度。
附图说明
图1为仿真计算模型的三维几何示意图；
图2为发射线框和接收测线的二维布局图；
图3为本实施例快速三维反演方法的流程图；
图4为仿真的瞬变电磁测量数据的CDI切面图；
图5为测量数据一阶矩、背景一阶矩、异常区域一阶矩在各条测线上曲线图；
图6为无约束条件的反演时间常数切面图；
图7为CDI起始模型的反演时间常数切面图；
图8为电导率加权值的切面图；
图9为采用电导率加权的反演时间常数切面图；
图10为深度加权值的切面图；
图11为采用深度加权的反演时间常数切面图；
图12为矩形大定源回线观测装置示意图。
具体实施方式
为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本发明进一步详细说明。需要说明的是，在附图或说明书描述中，相似或相同的部分都使用相同的图号。附图中未绘示或描述的实现方式，为所属技术领域中普通技术人员所知的形式。另外，虽然本文可提供包含特定值的参数的示范，但应了解，参数无需确切等于相应的值，而是可在可接受的误差容限或设计约束内近似于相应的值。实施例中提到的方向用语，例如“上”、“下”、“前”、“后”、“左”、“右”等，仅是参考附图的方向。因此，使用的方向用语是用来说明并非用来限制本发明的保护范围。
本发明采用了一种基于不同约束条件的瞬变电磁快速三维反演方法，通过简化的三维正演算法，将复杂的非线性严格反演问题转换成简单的线性反演问题；为了解决直接反演不能给出合理的结果，引入不同约束条件，快速、准确地实现了地下异常体的三维反演。
在本发明的一个示例性实施例中，提供了一种对带有Marco模块的Emit Maxwell软件的仿真瞬变电磁数据实施快速三维反演的方法。Marco模块由澳大利亚联邦科学工业研究院(CSIRO)开发，该模块基于三维积分方程，可以计算层状大地中含有多个棱柱体异常目标的瞬变电磁响应。
图1为仿真计算模型的三维几何示意图，如图所示，在电导率为1mS/m均匀大地背景中放置电导率为1S/m的平板异常体，平板大小为800E×800N×300Z(东西向800米×南北向800米×深度300米)，其中，E，N，Z分别代表东西向、南北向、深度向，该平板的上表面中心坐标为(0E，0N，-400Z)。
图3为本实施例快速三维反演方法的流程图。请参照图3，本实施例包括：
步骤A，在地面上铺设发射装置和接收点，设定地下异常区域，将异常区域划分为立方体形状的微元，根据发射装置、接收点以及微元的几何参数计算几何耦合因子矩阵；
步骤A具体包括：
子步骤A1：在被探测区域设置发射装置和接收点；
具体而言，该子步骤A1具体包括：在地面上铺设矩形大定源线框，线框中心位于(0E，0N，0Z)，线框边长为500米×500米，在南北向上从-500N到500N均匀分布11条接收测线，测线间距为100m，每条测线的走向从-1000E到1000E，均匀分布有21个接收点，共有接收点N＝231。
在本实施例中，发射装置为矩形大定源发射线框。图2为发射线框和接收测线的二维布局图。
子步骤A2：将异常区域划分为微元，其中，该异常区域为被探测区域内可能存在异常体的目标区域；具体而言，该子步骤A2具体包括：设定反演中异常区域范围为-1000E到1000E(东西向的-1000米至+1000米)，-500N到500N(南北方向的-500米至+500米)，-2000Z到-200Z(深度方向的-200米至-2000米)，将异常区域划分为K个边长为25m的立方体微元，微元总数量K＝230400；
子步骤A3：根据发射装置、接收点以及微元的几何参数计算几何耦合因子矩阵，几何耦合因子矩阵Gnk为：
上式中，表示第k个微元中心指向第n个接收点的单位方向矢量，表示入射到第k个微元的一次场的单位方向矢量，Vk表示第k个微元的体积，B0，k表示入射到第k个微元一次场幅度，rnk表示第k个微元到第n个接收点的距离。
本实施例中，子步骤A3的几何耦合因子矩阵Gnk大小为231×230400；
本实例中，几何耦合因子矩阵占用的储存空间为390MB。该几何耦合因子矩阵与瞬变电磁测量数据无关，它只与发射接收的几何位置信息有关，因此，该矩阵可以在测量之前就计算完成，在实施迭代反演中作为一个加载数据项；
步骤B：发射装置发射电流信号，电流关断后，各个接收点采集磁场数据，然后采用均匀大地模型，将采集的磁场数据转换为视电导率深度图(CDI)：
图4为仿真的瞬变电磁测量数据的CDI切面图，从图中可以看出，高电导率区域一定程度反应了平板异常体的位置，但是高电导率区域的底部位置远远超出了实际平板异常体的底部。从CDI中可以获取起始时间t1和截止时间tn处的视电导率σ1和σn，并估算背景电导率σbg。该CDI可以用于计算CDI起始模型和电导率加权值，在实地测量数据反演中，CDI也可以作为参考，验证反演结果的可靠性；
步骤C，根据瞬变电磁一阶矩变换，基于视电导率深度图(CDI)和测量磁场数据获取接收点处的异常区域的参考一阶矩；
瞬变电磁一阶矩变换定义为：
即磁场响应从0时刻到∞的积分。瞬变电磁阻性限制在时域中等于一阶矩，将磁场测量数据转换成阻性限制数据，等同一阶矩变换。
步骤C具体包括：
子步骤C1：
基于测量数据和均匀大地垂直磁场响应计算测量数据一阶矩；
子步骤C1具体包括：
由于测量数据在有限的时间范围内，为了得到测量数据的一阶矩，需要补齐时间范围外的磁场积分部分，这里，仿真的测量数据一阶矩的表达式为：
其中，t1和tn分别表示测量数据起始时间和截止时间，σ1和σn分别表示t1和tn处的视电导率，头部表示磁场从0到t1时间内的积分，中部表示测量数据从t1到tn的积分，尾部表示磁场从tn到∞时间内的积分。
中部是首先得到测量数据B(t)的拟合函数，将拟合函数在时间t1和tn内数值积分即得中部积分。由于在头部和尾部的时间范围内磁场无法测量，这里采用理论上均匀大地的磁场响应积分代替，头部通过计算从0到t1内均匀大地的磁场响应积分得到；尾部通过计算从tn到∞均匀大地的磁场响应积分得到，其中均匀大地的电导率从视电导率深度图中估计得到。
上述仿真的测量数据一阶矩中涉及理论上均匀大地的磁场响应积分的具体求解如下：
G(x，y，t)是单位电流元激发的均匀大地的垂直磁场响应，其表达式为：
其中，σ表示均匀大地的电导率，从视电导率深度图中估计得到，μ为真空磁导率，x和y表示相应接收点相对线电流元的末端的位置坐标，t为接收通道时间窗。
在本实施例中，由于采用矩形大定源回线，因此在矩形回线源激励下，均匀大地的磁场响应Bz(t)为：
上式中，G(x，y，t)是单位电流元的垂直磁场响应，x1＝XE-XR，x2＝XW-XR，y1＝YN-YR，y2＝YS-YR，其中(XE，YN)为矩形大定源回线东北方向的顶点坐标，(XW，YS)为矩形大定源回线西南方向的顶点坐标，如图12所示。
定义单位电流元的垂直磁场响应的积分形式：
当t→0时，上式中括号部分的表达式等于零，末端的积分部分可以进一步化简，因此，矩形大定源回线的其中一条边框激励形成的一阶矩的表达式为
上式表示从边框(x，y1)到(x，y2)的积分，在矩形大定源回线激励下，均匀大地的一阶矩等于四条边框的一阶矩之和。
矩形大定源回线的一条边框激励下，均匀大地的磁场响应从时间0到t1的积分为：
均匀大地的磁场响应从时间tn到∞的积分为：
本实施例中，根据表达式(8)计算矩形大定源回线每条边框的均匀大地的磁场响应从时间0到t1的积分，将四条边框的上述积分值求和，即可求得测量数据一阶矩中的头部；根据表达式(9)计算矩形大定源回线每条边框的均匀大地的磁场响应从时间tn到∞的积分，将四条边框的上述积分值求和，即可求得测量数据一阶矩中的尾部。
为了使仿真的测量数据贴近于实地测量数据，将上述仿真的测量数据一阶矩加入5％的高斯白噪声；
子步骤C2：从视电导率深度图(CDI)中估计背景电导率，然后利用均匀大地垂直磁场的一阶矩表达式求出背景响应一阶矩；
子步骤C2具体包括：从视电导率深度图CDI中获取起始时间t1和截止时间tn处的视电导率σ1和σn，并估算背景电导率σbg。
其中上述估算背景电导率σbg具体包括：
背景响应一阶矩可由均匀大地的磁场响应从0到∞时间内积分得到。
在本实施例中，背景响应一阶矩根据表达式(7)，将矩形大定源回线的四条边框的一阶矩相加得到。
子步骤C3：从测量数据一阶矩中剔除背景一阶矩得到异常区域的参考一阶矩；
图5为测量数据一阶矩、背景一阶矩、异常区域一阶矩在各条测线上曲线图。该异常区域的参考一阶矩作为接下来的反演中的参考一阶矩d＝(d1，d2，…，dN)T，该实例中，接收点数N＝231。
步骤D，构建反演中时间常数向量的约束条件；所述步骤D中的约束条件为以下三种约束条件中的一种：
约束条件a：在反演过程中，时间常数向量中的值限制为正值，在反演起始时，时间常数向量的起始值由视电导率深度图得到；
约束条件b：在反演过程中，时间常数向量中的值限制为正值，并根据微元的视电导率给时间常数向量中的值赋予加权值，在反演起始时，时间常数向量的起始值为零；
约束条件c：在反演过程中，时间常数向量中的值限制为正值，并根据微元的深度给时间常数向量中的值赋予加权值，反演起始时，时间常数向量的起始值为零；
其中，
所述约束条件a中的时间常数向量的起始值由以下方式获得：
首先，定义起始时间常数向量第k个微元的时间常数为
其中，σk为第k个微元的视电导率值，由视电导率深度图CDI获得，μ为真空磁导率，L为微元的边长。
根据最小二乘法，实测响应与起始理论响应的拟合差为
其中β为拟合系数，dn为异常区域的参考一阶矩。为使拟合差最小，有偏微分方程求解该偏微分方程，得到β，起始时间常数向量为τ＝β·τ0。
所述约束条件b中加权值计算方法为：
电导率加权值与微元在CDI中的位置相关，其分布于0和1之间，第k个微元电导率加权值为
其中，σk为第k个微元的视电导率值，由视电导率深度图CDI获得，σmax为整个CDI中最大的视电导率。电导率加权使CDI中大电导率区域占主导作用，减弱小电导率区域的几何耦合因子影响。
所述约束条件c的加权值计算方法为：
深度加权值与微元的深度相关，第k个微元的深度加权值为
W(zk，s0，d0)＝tanh(s0·(zk-d0)) (13)
上式中，zk为第k个微元的深度，s0为斜率因子，d0为参考深度。对不同深度微元的时间常数赋予加权值，深度加权减弱了浅层耦合因子的影响。
步骤E，基于异常区域的参考一阶矩及其误差、微元对应的时间常数向量以及几何耦合因子矩阵构建反演运算的目标函数；
步骤E具体包括：
cn表示第n接收点处异常区域正演运算的理论一阶矩，其计算方法为：
理论上，大地总一阶矩响应由背景响应一阶矩和异常体响应一阶矩叠加构成，表达式如下：
上式中，表示第n个接收点处的背景响应一阶矩，等式右边求和部分表示异常体响应一阶矩，Gnk为第k个微元相对于第n个接收点的几何耦合因子矩阵，τk为第k个微元的时间常数。大地的总一阶矩理论上应等于测量数据一阶矩，在这里，需从测量数据一阶矩中反演出未知的τk。
为简化运算，在反演过程中背景响应一阶矩已从测量数据一阶矩中剔除，异常区域正演运算的理论一阶矩取：
则反演运算的目标函数表达式为
其中，微元的时间常数向量为τ＝(τ1，τ2，…，τK)T，τk为第k个微元的时间常数，异常区域的参考一阶矩的误差为q＝(q1，q2，…，qN)T，qn为第n个异常区域的参考一阶矩的误差，cn表示第n接收点处异常区域正演运算的理论一阶矩，Gnk为几何耦合因子矩阵，dn为异常区域的参考一阶矩。
步骤F，基于微元对应的时间常数向量的约束条件，采用最优化算法对反演运算的目标函数进行迭代；
步骤F具体包括：
子步骤F1：
若采用约束条件a，则第i+1次迭代前的时间常数向量为τi，第i+1次迭代后的时间常数向量为τi+1＝τi+δτi；
若采用约束条件b，则第i+1次迭代前的时间常数向量为τi，第i+1次迭代后的时间常数向量为τi+1＝τi+W(σk)δτi；
若采用约束条件c，则第i+1次迭代前的时间常数向量为τi，第i+1次迭代后的时间常数向量为τi+1＝τi+W(zk，s0，d0)δτi；
子步骤F2：
采用最优化算法对所述目标函数进行迭代。
最优化算法为以下方法中的一种：最速下降法、共轭梯度法、牛顿法、拟牛顿法、信赖域法、奇异值分解法、模拟退火法、遗传算法。
本实施例中，以最速下降法为例说明迭代方法，
δτi通过下式计算：
步长αi的表达式为
上式中，中间变量和的表达式为
其中，是拟合差目标函数的梯度，表示为
步骤G，判断迭代后的目标函数是否收敛，若收敛，则保存最优化时间常数向量作为反演最终结果，即τ＝(τ1，τ2，…，τK)T，否则，返回步骤F执行，其中，τk是异常区域内第k个微元的电导率相对应的时间常数，根据时间常数向量的反演最终结果得到异常体的位置。
其中，所述目标函数是否收敛是指目标函数是否小于1，小于1则收敛，否则不收敛。
时间常数τk与相应的第k个微元的电导率有关，电导率越大，时间常数越大。
图6为无约束条件的反演时间常数切面图，该反演结果经过73次迭代，耗时60s。由于反演过程中无约束条件，反演时间常数出现负值，同时，因为浅层的几何耦合因子更大，反演结果中异常区域贴近地面。可见无约束反演并不能反映平板异常体的真实地电结构，因此，反演过程中，设定约束条件是非常必要的。
图7为步骤D中的约束条件a涉及的CDI起始模型的反演时间常数切面图。
CDI起始模型的初始拟合差为1.58，迭代8次，时间0.5s后达到拟合差小于为1的收敛标准。从图7可以看出，时间常数反演结果会趋近CDI形态，反演的时间常数最大值与平板异常体的真实时间常数值(τslab＝4.5msec)非常吻合。
图8为步骤D中的约束条件b涉及的电导率加权值的切面图。
电导率加权使CDI中大电导率区域占主导作用，减弱小电导率区域的几何耦合因子影响。
图9为采用电导率加权的反演时间常数切面图。初始拟合差为67，反演迭代50次后达到收敛，耗时43s，如图所示，不同于CDI起始模型的反演结果，电导率加权反演结果中高时间常数区域局限于实际平板的体积范围内，较准确地反映了异常体的位置和体积。
图10为步骤D中的约束条件c涉及的深度加权值的切面图。
本次实例中，s0＝0.004，d0＝350，该设定参数计算出的深度加权值如图10所示。
图11为采用深度加权的反演时间常数切面图。初始拟合差为67，迭代5次后达到收敛，耗时8s，如图所示，高时间常数区域与平板异常体的位置和体积较好的吻合。
至此，已经结合附图对本实施例进行了详细描述。依据以上描述，本领域技术人员应当对本发明基于不同约束条件的瞬变电磁快速三维反演方法有了清楚的认识。
此外，上述对各元件和方法的定义并不仅限于实施例中提到的各种具体结构、形状或方式，本领域普通技术人员可对其进行简单地更改或替换。
综上所述，本发明提供了一种基于不同约束条件的瞬变电磁快速三维反演方法，该方法采用了一种数据压缩的处理方法及简化的三维正演算法，解决了当前三维反演面临的数据量大、正演复杂的难题；由于三维反演问题的欠定性，直接反演一般不能给出合理的结果，为了促使反演结果贴近实际地下结构，本方法提出了不同的约束条件，可以较为准确反演出地下异常体。本发明提出的方法可在普通计算机上实现地下异常目标的快速三维反演，在实时瞬变电磁数据解释工作中具有重要应用前景。
以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。
Claims (9)
1.一种基于不同约束条件的瞬变电磁快速三维反演方法，其特征在于，包括：
步骤A：在地面上铺设发射装置和接收点，设定地下异常区域，将异常区域划分为立方体形状的微元，根据发射装置、接收点以及微元的几何参数计算几何耦合因子矩阵；所述几何耦合因子矩阵Gnk为：
<mrow>
<msub>
<mi>G</mi>
<mrow>
<mi>n</mi>
<mi>k</mi>
</mrow>
</msub>
<mo>=</mo>
<msub>
<mi>B</mi>
<mrow>
<mn>0</mn>
<mo>,</mo>
<mi>k</mi>
</mrow>
</msub>
<msub>
<mi>V</mi>
<mi>k</mi>
</msub>
<mfrac>
<mrow>
<mn>3</mn>
<mrow>
<mo>(</mo>
<msub>
<mover>
<mi>b</mi>
<mo>^</mo>
</mover>
<mi>k</mi>
</msub>
<mo>&times;</mo>
<msub>
<mover>
<mi>r</mi>
<mo>^</mo>
</mover>
<mrow>
<mi>n</mi>
<mi>k</mi>
</mrow>
</msub>
<mo>)</mo>
</mrow>
<msub>
<mover>
<mi>r</mi>
<mo>^</mo>
</mover>
<mrow>
<mi>n</mi>
<mi>k</mi>
</mrow>
</msub>
<mo>-</mo>
<msub>
<mover>
<mi>b</mi>
<mo>^</mo>
</mover>
<mi>k</mi>
</msub>
</mrow>
<mrow>
<mn>4</mn>
<msup>
<msub>
<mi>&pi;r</mi>
<mrow>
<mi>n</mi>
<mi>k</mi>
</mrow>
</msub>
<mn>3</mn>
</msup>
</mrow>
</mfrac>
<mfrac>
<msup>
<mi>&pi;</mi>
<mn>2</mn>
</msup>
<mn>10</mn>
</mfrac>
</mrow>
上式中，表示第k个微元中心指向第n个接收点的单位方向矢量，表示入射到第k个微元的一次场的单位方向矢量，Vk表示第k个微元的体积，B0,k表示入射到第k个微元一次场幅度，rnk表示第k个微元到第n个接收点的距离；
步骤B：发射装置发射电流信号，电流关断后，各个接收点采集磁场数据，然后采用均匀大地模型，将采集的磁场数据转换为视电导率深度图；
步骤C：根据瞬变电磁一阶矩变换，基于视电导率深度图和测量磁场数据获取接收点处的异常区域的参考一阶矩；
步骤D：构建反演中时间常数向量的约束条件；
步骤E：基于异常区域的参考一阶矩及其误差、微元对应的时间常数向量以及几何耦合因子矩阵构建反演运算的目标函数；
步骤F：基于微元对应的时间常数向量的约束条件，采用最优化算法对反演运算的目标函数进行迭代；以及
步骤G：判断迭代后的目标函数是否收敛，若收敛，则保存最优化时间常数向量作为反演最终结果，即τ＝(τ1,τ2,…,τK)T，否则，返回步骤F执行，其中，τk是异常区域内第k个微元的电导率对应的时间常数，根据时间常数向量的反演最终结果得到异常体的位置和体积。
2.根据权利要求1所述的瞬变电磁快速三维反演方法，其特征在于，所述步骤D中的约束条件为以下三种约束条件中的一种：
约束条件a：在反演过程中，时间常数向量中的值限制为正值，在反演起始时，时间常数向量的起始值由视电导率深度图得到；
约束条件b：在反演过程中，时间常数向量中的值限制为正值，并根据微元的视电导率给时间常数向量中的值赋予加权值，在反演起始时，时间常数向量的起始值为零；
约束条件c：在反演过程中，时间常数向量中的值限制为正值，并根据微元的深度给时间常数向量中的值赋予加权值，反演起始时，时间常数向量的起始值为零。
3.根据权利要求2所述的瞬变电磁快速三维反演方法，其特征在于，所述约束条件a中的时间常数向量的起始值由以下方式获得
定义起始时间常数向量第k个微元的时间常数为：
<mrow>
<msubsup>
<mi>&tau;</mi>
<mi>k</mi>
<mn>0</mn>
</msubsup>
<mo>=</mo>
<mfrac>
<mrow>
<msub>
<mi>&sigma;</mi>
<mi>k</mi>
</msub>
<msup>
<mi>&mu;L</mi>
<mn>2</mn>
</msup>
</mrow>
<mn>25.6</mn>
</mfrac>
</mrow>
其中，σk为第k个微元的视电导率值，由视电导率深度图获得，μ为真空磁导率，L为微元的边长；
根据最小二乘法，实测响应与起始理论响应的拟合差为：
<mrow>
<mi>S</mi>
<mo>=</mo>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>n</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>N</mi>
</munderover>
<msup>
<mrow>
<mo>(</mo>
<msub>
<mi>d</mi>
<mi>n</mi>
</msub>
<mo>-</mo>
<mi>&beta;</mi>
<mo>&CenterDot;</mo>
<msubsup>
<mi>c</mi>
<mi>n</mi>
<mn>0</mn>
</msubsup>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
</mrow>
其中β为拟合系数，dn为异常区域的参考一阶矩，Gnk为几何耦合因子矩阵，为使拟合差最小，有偏微分方程求解该偏微分方程，得到β，起始时间常数向量为τ＝β·τ0。
4.根据权利要求2所述的瞬变电磁快速三维反演方法，其特征在于，所述约束条件b中加权值计算方法为：
电导率加权值与微元在视电导率深度图中的位置相关，其分布于0和1之间，第k个微元电导率加权值为：
<mrow>
<mi>W</mi>
<mrow>
<mo>(</mo>
<msub>
<mi>&sigma;</mi>
<mi>k</mi>
</msub>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mfrac>
<msub>
<mi>&sigma;</mi>
<mi>k</mi>
</msub>
<msub>
<mi>&sigma;</mi>
<mrow>
<mi>m</mi>
<mi>a</mi>
<mi>x</mi>
</mrow>
</msub>
</mfrac>
</mrow>
其中，σk为第k个微元的视电导率值，由视电导率深度图获得，σmax为整个视电导率深度图中最大的视电导率。
5.根据权利要求2所述的瞬变电磁快速三维反演方法，其特征在于，所述约束条件c的加权值计算方法为：
深度加权值与微元的深度相关，第k个微元的深度加权值为：
W(zk,s0,d0)＝tanh(s0·(zk-d0))
上式中，zk为第k个微元的深度，s0为斜率因子，d0为参考深度。
6.根据权利要求2所述的瞬变电磁快速三维反演方法，其特征在于，所述步骤F具体包括：
子步骤F1：
若采用约束条件a，则第i+1次迭代前的时间常数向量为τi，第i+1次迭代后的时间常数向量为τi+1＝τi+δτi；
若采用约束条件b，则第i+1次迭代前的时间常数向量为τi，第i+1次迭代后的时间常数向量为τi+1＝τi+W(σk)δτi；
若采用约束条件c，则第i+1次迭代前的时间常数向量为τi，第i+1次迭代后的时间常数向量为τi+1＝τi+W(zk,s0,d0)δτi；以及
子步骤F2：
采用最优化算法对所述目标函数进行迭代。
7.根据权利要求1至6中任一项权利要求所述的瞬变电磁快速三维反演方法，其特征在于，所述步骤A具体包括：
子步骤A1：在被探测区域设置发射装置和接收点；
子步骤A2：将异常区域划分为微元，其中，该异常区域为被探测区域内可能存在异常体的目标区域；以及
子步骤A3：根据发射装置、接收点以及微元的几何参数计算几何耦合因子矩阵。
8.根据权利要求7所述的瞬变电磁快速三维反演方法，其特征在于，所述步骤E具体包括：
反演运算的目标函数表达式为：
<mrow>
<mi>S</mi>
<mrow>
<mo>(</mo>
<mi>&tau;</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mfrac>
<mn>1</mn>
<mi>N</mi>
</mfrac>
<munderover>
<mi>&Sigma;</mi>
<mrow>
<mi>n</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>N</mi>
</munderover>
<msup>
<mrow>
<mo>(</mo>
<mfrac>
<mrow>
<msub>
<mi>d</mi>
<mi>n</mi>
</msub>
<mo>-</mo>
<msub>
<mi>c</mi>
<mi>n</mi>
</msub>
</mrow>
<msub>
<mi>q</mi>
<mi>n</mi>
</msub>
</mfrac>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
</mrow>
其中，微元的时间常数向量为τ＝(τ1,τ2,…,τK)T，τk为第k个微元的时间常数，异常区域的参考一阶矩的误差为q＝(q1,q2,…,qN)T，qn为第n个异常区域的参考一阶矩的误差，cn表示第n接收点处异常区域正演运算的理论一阶矩，Gnk为几何耦合因子矩阵，dn为异常区域的参考一阶矩。
9.根据权利要求1至6中任一项权利要求所述的瞬变电磁快速三维反演方法，其特征在于，所述步骤C具体包括：
子步骤C1：基于测量数据和均匀大地垂直磁场响应计算测量数据一阶矩；
子步骤C2：从视电导率深度图中估计背景电导率，然后利用均匀大地垂直磁场的一阶矩表达式求出背景响应一阶矩；以及
子步骤C3：从测量数据一阶矩中剔除背景一阶矩得到异常区域的参考一阶矩。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201510926452.7A CN105589108B (zh) | 2015-12-14 | 2015-12-14 | 基于不同约束条件的瞬变电磁快速三维反演方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201510926452.7A CN105589108B (zh) | 2015-12-14 | 2015-12-14 | 基于不同约束条件的瞬变电磁快速三维反演方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN105589108A CN105589108A (zh) | 2016-05-18 |
| CN105589108B true CN105589108B (zh) | 2017-11-21 |
Family
ID=55928832
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201510926452.7A Active CN105589108B (zh) | 2015-12-14 | 2015-12-14 | 基于不同约束条件的瞬变电磁快速三维反演方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN105589108B (zh) |
Families Citing this family (13)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106646625A (zh) * | 2016-09-27 | 2017-05-10 | 中国科学院电子学研究所 | 一种锐边界模型的瞬变电磁反演方法 |
| CN108802834B (zh) * | 2018-02-13 | 2020-12-22 | 中国科学院电子学研究所 | 一种基于联合反演的地下目标识别方法 |
| CN110222429A (zh) * | 2019-06-10 | 2019-09-10 | 清华大学 | 多根线电流参数重建的优化方法 |
| CN110865414B (zh) * | 2019-11-01 | 2021-02-02 | 吉林大学 | 一种用于城市地下空间探测的瞬变电磁噪声抑制方法 |
| CN110989002B (zh) * | 2019-11-07 | 2021-01-05 | 吉林大学 | 地空时域电磁***浅部低阻异常体数据解释方法 |
| CN111796335B (zh) * | 2020-08-28 | 2023-03-14 | 核工业航测遥感中心 | 航空瞬变电磁时间常数提取方法 |
| CN112965120B (zh) * | 2021-02-01 | 2024-04-12 | 西双版纳景海高速公路建设投资有限公司 | 一种地质勘探用瞬变电磁信号处理方法、装置及存储介质 |
| CN112949134B (zh) * | 2021-03-09 | 2022-06-14 | 吉林大学 | 基于非结构有限元方法的地-井瞬变电磁反演方法 |
| CN113176617A (zh) * | 2021-03-15 | 2021-07-27 | 中煤科工集团西安研究院有限公司 | 一种沉积地层瞬变电磁多参数约束反演成像方法 |
| CN113552637B (zh) * | 2021-07-30 | 2023-11-14 | 中国自然资源航空物探遥感中心 | 一种航空-地面-井中磁异常数据协同三维反演方法 |
| CN113391362B (zh) * | 2021-08-13 | 2021-10-29 | 成都理工大学 | 基于廊带数据约束的大地电磁剖面三维结构化反演方法 |
| CN114114438B (zh) * | 2021-09-27 | 2023-07-18 | 中国地质科学院地球物理地球化学勘查研究所 | 一种回线源地空瞬变电磁数据的拟三维反演方法 |
| CN114265124B (zh) * | 2021-12-30 | 2022-08-02 | 成都理工大学 | 基于时间域瞬变电磁概率反演的不良地质体定位方法 |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103760614A (zh) * | 2014-02-24 | 2014-04-30 | 中国科学院电子学研究所 | 一种适用于不规则发射波形的瞬变电磁正演方法 |
| CN103777248A (zh) * | 2014-02-08 | 2014-05-07 | 中国科学院电子学研究所 | 一种适用于不规则发射回线的tem一维正演方法 |
| CN104360404A (zh) * | 2014-11-27 | 2015-02-18 | 中国科学院电子学研究所 | 基于不同约束条件的大地电磁正则化反演方法 |
Family Cites Families (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US8095345B2 (en) * | 2009-01-20 | 2012-01-10 | Chevron U.S.A. Inc | Stochastic inversion of geophysical data for estimating earth model parameters |
| CN102419454A (zh) * | 2011-06-30 | 2012-04-18 | 中国科学院地质与地球物理研究所 | 隧道掌子面前方远距离含水目标体的瞬变电磁预报方法 |
| CN102901989A (zh) * | 2011-07-29 | 2013-01-30 | 中国国土资源航空物探遥感中心 | 一种基于重力场或磁场数据的地质体三维可视化建模与解释方法 |
-
2015
- 2015-12-14 CN CN201510926452.7A patent/CN105589108B/zh active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103777248A (zh) * | 2014-02-08 | 2014-05-07 | 中国科学院电子学研究所 | 一种适用于不规则发射回线的tem一维正演方法 |
| CN103760614A (zh) * | 2014-02-24 | 2014-04-30 | 中国科学院电子学研究所 | 一种适用于不规则发射波形的瞬变电磁正演方法 |
| CN104360404A (zh) * | 2014-11-27 | 2015-02-18 | 中国科学院电子学研究所 | 基于不同约束条件的大地电磁正则化反演方法 |
Non-Patent Citations (1)
| Title |
|---|
| 时间域电磁勘探数据的模拟退火法反演研究;张建国，等;《电子与信息学报》;20150131;第37卷(第1期);220-225 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN105589108A (zh) | 2016-05-18 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN105589108B (zh) | 基于不同约束条件的瞬变电磁快速三维反演方法 | |
| Grayver et al. | 3D inversion and resolution analysis of land-based CSEM data from the Ketzin CO 2 storage formation | |
| Yang et al. | Three-dimensional inversion of airborne time-domain electromagnetic data with applications to a porphyry deposit | |
| Hansen et al. | Accounting for imperfect forward modeling in geophysical inverse problems—exemplified for crosshole tomography | |
| CN101609169B (zh) | 一种提高电磁波电阻率测量精度和扩展其测量范围的方法 | |
| CN112949134B (zh) | 基于非结构有限元方法的地-井瞬变电磁反演方法 | |
| CN104360401B (zh) | 一种瞬变电磁b场确定地下目标体地质信息方法 | |
| CN106291725B (zh) | 一种快速反演地下地质体空间位置的方法 | |
| CN105386756B (zh) | 一种应用应变量计算脆性地层孔隙度的方法 | |
| CN106199742A (zh) | 一种频率域航空电磁法2.5维带地形反演方法 | |
| WO2009146041A1 (en) | Constructing a reduced order model of an electromagnetic response in a subterranean structure | |
| CN104854480A (zh) | 用以在地下岩层中寻找位置的设备和方法 | |
| CN103728667A (zh) | 一种视三维高密度电法的浅表层地质结构建模方法 | |
| CN103777248A (zh) | 一种适用于不规则发射回线的tem一维正演方法 | |
| CN110458129A (zh) | 基于深度卷积神经网络的非金属地雷识别方法 | |
| CN105573963B (zh) | 一种电离层水平不均匀结构重构方法 | |
| CN104267443B (zh) | 基于反演模型的大地电磁场静位移校正方法 | |
| CN104316961A (zh) | 获取风化层的地质参数的方法 | |
| Zhang et al. | 3D inversion of large-scale frequency-domain airborne electromagnetic data using unstructured local mesh | |
| Geng et al. | 3D joint inversion of gravity-gradient and borehole gravity data | |
| CN116879964A (zh) | 一种时频电磁频率域数据自约束稳健电阻率反演方法 | |
| CN104020508A (zh) | 一种用于地质雷达层析探测的直射线追踪算法 | |
| CN105550442B (zh) | 基于瞬变电磁矩变换的数据处理及三维正演方法 | |
| CN114114438B (zh) | 一种回线源地空瞬变电磁数据的拟三维反演方法 | |
| CN115469365A (zh) | 一种堤防渗漏探测方法、装置、电子设备以及存储介质 |
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 |