CN114970290A - 一种地面瞬变电磁法并行反演方法及系统 - Google Patents
一种地面瞬变电磁法并行反演方法及系统 Download PDFInfo
- Publication number
- CN114970290A CN114970290A CN202210888762.4A CN202210888762A CN114970290A CN 114970290 A CN114970290 A CN 114970290A CN 202210888762 A CN202210888762 A CN 202210888762A CN 114970290 A CN114970290 A CN 114970290A
- Authority
- CN
- China
- Prior art keywords
- model
- transient electromagnetic
- inversion
- response data
- electromagnetic response
- 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.)
- Granted
Links
- 230000001052 transient effect Effects 0.000 title claims abstract description 193
- 238000000034 method Methods 0.000 title claims abstract description 165
- 230000004044 response Effects 0.000 claims abstract description 114
- 239000011159 matrix material Substances 0.000 claims abstract description 71
- 230000035945 sensitivity Effects 0.000 claims abstract description 64
- 239000013598 vector Substances 0.000 claims abstract description 44
- 238000011156 evaluation Methods 0.000 claims abstract description 9
- 238000010606 normalization Methods 0.000 claims abstract description 7
- 208000037516 chromosome inversion disease Diseases 0.000 claims description 174
- 230000005684 electric field Effects 0.000 claims description 65
- 238000004458 analytical method Methods 0.000 claims description 38
- 238000004364 calculation method Methods 0.000 claims description 28
- 238000004088 simulation Methods 0.000 claims description 7
- 238000010276 construction Methods 0.000 claims description 5
- 238000009795 derivation Methods 0.000 claims description 5
- 238000013461 design Methods 0.000 claims description 5
- 230000017105 transposition Effects 0.000 claims 1
- 230000006870 function Effects 0.000 description 17
- 230000008569 process Effects 0.000 description 13
- 238000000354 decomposition reaction Methods 0.000 description 10
- 238000005259 measurement Methods 0.000 description 9
- 238000004422 calculation algorithm Methods 0.000 description 6
- 238000005516 engineering process Methods 0.000 description 6
- 230000005540 biological transmission Effects 0.000 description 5
- 238000010586 diagram Methods 0.000 description 5
- 230000002159 abnormal effect Effects 0.000 description 4
- 238000001514 detection method Methods 0.000 description 4
- 238000006467 substitution reaction Methods 0.000 description 3
- 238000012937 correction Methods 0.000 description 2
- 230000007547 defect Effects 0.000 description 2
- 238000003384 imaging method Methods 0.000 description 2
- 230000003044 adaptive effect Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 238000001816 cooling Methods 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 238000009792 diffusion process Methods 0.000 description 1
- 238000006073 displacement reaction Methods 0.000 description 1
- 238000005553 drilling Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 230000005672 electromagnetic field Effects 0.000 description 1
- 230000007613 environmental effect Effects 0.000 description 1
- 238000013210 evaluation model Methods 0.000 description 1
- 230000005284 excitation Effects 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 230000035772 mutation Effects 0.000 description 1
- BULVZWIRKLYCBC-UHFFFAOYSA-N phorate Chemical compound CCOP(=S)(OCC)SCSCC BULVZWIRKLYCBC-UHFFFAOYSA-N 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000005070 sampling Methods 0.000 description 1
- 238000010183 spectrum analysis Methods 0.000 description 1
- 238000004613 tight binding model Methods 0.000 description 1
Images
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
- G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
-
- 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/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
- G06F17/12—Simultaneous equations, e.g. systems of linear equations
-
- 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/30—Circuit design
- G06F30/36—Circuit design at the analogue level
- G06F30/367—Design verification, e.g. using simulation, simulation program with integrated circuit emphasis [SPICE], direct methods or relaxation methods
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
- G06T17/20—Finite element generation, e.g. wire-frame surface description, tesselation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/12—Timing analysis or timing optimisation
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
- Y02A90/30—Assessment of water resources
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Computer Hardware Design (AREA)
- Mathematical Physics (AREA)
- General Engineering & Computer Science (AREA)
- Geometry (AREA)
- Mathematical Optimization (AREA)
- Mathematical Analysis (AREA)
- Computational Mathematics (AREA)
- Pure & Applied Mathematics (AREA)
- Evolutionary Computation (AREA)
- Data Mining & Analysis (AREA)
- Software Systems (AREA)
- Algebra (AREA)
- Databases & Information Systems (AREA)
- Operations Research (AREA)
- Computing Systems (AREA)
- Computer Graphics (AREA)
- Microelectronics & Electronic Packaging (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
本发明提供了一种地面瞬变电磁法并行反演方法及系统，属于地球物理方法技术领域，采集观测瞬变电磁响应数据；将反演空间剖分成非结构化四面体网格并设置电导率值，构建初始反演模型以及正则化目标函数；计算模型参数的灵敏度矩阵与向量的乘积；采用高速‑牛顿法方程将目标函数转化为最小二乘问题，获取模型更新方向，并采用线性搜索获得最优模型更新步长，更新反演模型；基于矢量有限元方法进行瞬变电磁三维正演，获取预测瞬变电磁响应数据；采用归一化误差评价预测瞬变电磁响应数据与观测瞬变电磁响应数据的拟合情况；若数据归一化拟合差达到预设阈值反演方法终止。本发明可以精确模拟复杂地质结构的瞬变电磁响应。
Description
技术领域
本发明属于地球物理方法技术领域，更具体地，涉及一种地面瞬变电磁法并行反演方法及系统。
背景技术
电导率是地球内部的物性特征之一，获取地下电阻率的分布信息可以为解决地球科学问题和工程问题提供关键支撑。瞬变电磁法作为地球物理电磁法的主要分支，是获取地下电阻率分布信息的重要手段。瞬变电磁法具有频带宽、抗干扰能力强和探测装置灵活等优点，在能源资源、环境工程、地灾监测等领域取得广泛应用。因此，瞬变电磁数据精细可靠的反演地下电阻率模型极为重要。
目前，瞬变电磁数据解释通常采用视电阻率成像和一维反演方法进行瞬变电磁数据解释。前者可用于定性解释，但视电阻率的定义是模糊不精确的，该方法只适用于某些地质条件简单的地区。将地下结构简化为一维的层状模型，采用一维反演可以获得层状电导率模型，该方法是瞬变电磁数据解释的主流方法之一，但根据一维反演结果构建的拟二维剖面通常会产生电阻率突变以至于难以进行地质解释。在一维反演中引入横向约束和空间约束技术，一维反演结果的空间连续性一定程度上得到了增强。但在复杂的地质情况下，视电阻率成像方法和一维反演不能恢复精确的地下结构。
三维反演是对瞬变电磁数据精确解释的有效方法。正演是反演的基础，积分方程法和有限差分法被应用于电磁三维正演模拟。不同的方法优缺点不同。有限差分法原理简单，容易算法实现，但难以对复杂几何结构进行精确剖分是其难以弥补的缺点。积分方程法的模拟精度严重依赖系数矩阵和并矢格林函数解的精度，对于大规模复杂问题求解线性方程组对内存的需求极大增加，求解速度慢且求解精度降低。有限元方法可以通过使用灵活的非结构化网格剖分技术模拟复杂地形、地质体以及任意形态的反射源，具有非常高的计算精度，最适合应用于精细探测，但计算效率仍有待提高。有限体积法也可以采用非结构网格离散研究区域，但其精度仍然比不上有限元方法。
基于规则网格的积分方程法、有限差分法、有限体积法和瞬变电磁三维反演也有少量报道，但采用了规则网格对模型进行离散，难以考虑复杂的地质构造（如地形）的影响。同时，瞬变电磁法三维反演方法的计算需求非常庞大，在很大程度上限制了瞬变电磁数据三维反演在实际中的应用。
综上所述，为了实现瞬变电磁数据的精确解释，有必要解决瞬变电磁三维正反演的地形影响和计算效率问题，使瞬变电磁精细探测成为可能。
发明内容
针对现有技术的缺陷，本发明的目的在于提供一种地面瞬变电磁法并行反演方法及系统，旨在解决现有的瞬变电磁三维正反演缺乏考虑地形影响且计算效率问题的问题。
为实现上述目的，一方面，本发明提供了一种地面瞬变电磁法并行反演方法，包括以下步骤：
S1：获取用于反演的观测瞬变电磁响应数据，包括采集地面的实际观测瞬变电磁响应数据和采用理论模型通过瞬变电磁三维正演获取理论观测瞬变电磁响应数据；理论模型为能够获取理论观测瞬变电磁响应数据的模拟模型；
S2：将反演空间剖分成非结构化四面体网格；其中，反演空间为测点及其附近预设延伸范围所对应的地下空间；
S3：在非结构化四面体网格中设置电导率值，构建初始反演模型，并结合观测瞬变电磁响应数据构建初始反演模型的正则化目标函数；
S4：对棱边电场方程中初始反演模型中的模型参数进行求导，采用反向时间步进的隐式方法，逐步计算电场沿所有单元棱边相对模型参数的灵敏度转置，再结合系数插值矩阵，计算观测瞬变电磁响应数据对模型参数的灵敏度矩阵与向量的乘积；
S5：结合观测瞬变电磁响应数据对模型参数的灵敏度与向量的乘积，采用具有拟二次收敛速度的高速-牛顿法方程将目标函数转化为最小二乘问题，获取模型更新方向，并采用线性搜索获得最优模型更新步长；
S6：结合模型更新方向与最优模型更新步长计算模型参数，更新反演模型；
S7：基于矢量有限元方法对更新后的反演模型进行瞬变电磁三维正演，获取预测瞬变电磁响应数据；
S8：采用归一化误差评价预测瞬变电磁响应数据与观测瞬变电磁响应数据的拟合情况；若数据归一化拟合差达到预设阈值，反演方法终止；若未达到预设阈值，且当前迭代次数小于设定次数，则转至S4。
进一步优选地，获取理论观测瞬变电磁响应数据的方法，包括以下步骤：
对时域的准静态Maxwell方程组两边取旋度，并采用基于非结构化四面体网格的一阶矢量有限元方法，获取有限元分析方程；其中，采用网格单元的棱边离散瞬变电磁发射源，使导线或线圈与网格单元的棱边重合；
采用后退二阶欧拉方法对有限元分析方程中的时间进行离散，且通过泰勒级数展开，获取各时刻非均匀步长的电场时间导数表达式，对有限元分析方程进行变形；
采用Dirichlet边界条件，假设电场在建模区域的边界上消失，设置电场为零的初始条件，且将理论模型的建模区域扩展，在测点和发射源的预设区域加密网格，其他区域逐步增大网格的尺寸，完成边界条件和初始条件的设定以及网格的设置；
基于边界条件和初始条件的设定以及网格的设置，采用直接求解器求解变形后有限元分析方程，计算各时刻的电场；
根据电场和稀疏插值矩阵计算理论观测瞬变电磁响应数据。
进一步优选地，正则化目标函数为：
φ（m）=φd（m）+λφm（m）
φd（m）=0.5||Wd（F（m）-dobs）||2 L2
φm（m）=0.5||Wm（m-mref）||2 L2
其中，φd是数据拟合项；φm是模型稳定项；λ是正则化参数；m为反演模型参数；m=log（σ）；σ为电导率；Wd=diag（1/（dobs+ε））；F（m）表示初始反演模型参数m的正演响应， dobs为测点上的观测瞬变电磁响应数据，Wm为模型粗糙度矩阵，mref为含先验信息的参考模型参数；x是随机向量；J为观测瞬变电磁响应数据对模型参数的灵敏度；L2为二范数；ε为预设的较小值。
进一步优选地，Sn=∂En/∂m；
其中，S n 为第n个时间步电场沿所有单元棱边相对于模型参数的灵敏度；S n 的大小是N ed ×N m ，N ed 是网格棱边的数量；N m 是模型参数的数量；
Jn=QSn；
其中，Q是大小为N d×N ed的稀疏插值矩阵；Jn为第n个时间步观测瞬变电磁响应数据对模型参数的灵敏度；N d表示数据个数；数据个数为测取测点数据所用仪器对应时间道的数目与测点数目的乘积。
进一步优选地，m k =m k-1+lδm k
其中，m k 为第k次反演模型的模型参数；m k-1为第k-1次反演模型的模型参数；δmk为模型更新方向；l为模型更新步长。
另一方面，本发明提供了一种地面瞬变电磁法并行反演系统，其特征在于，包括：
瞬变电磁数据的获取模块，用于获得用于反演的观测瞬变电磁响应数据，包括野外实地数据采集单元与理论模型瞬变电磁三维正演单元，分别获取实际观测瞬变电磁响应数据和理论观测瞬变电磁响应数据；
空间剖分模块，用于将反演空间剖分成非结构化四面体网格；其中，反演空间为测点附近的预设地下空间；
初始反演模型的构建模块，用于在非结构化四面体网格中设置电导率值，构建初始反演模型，并结合观测瞬变电磁响应数据构建初始反演模型的正则化目标函数；
观测瞬变电磁响应数据对模型参数的灵敏度计算模块，用于对棱边电场方程中初始反演模型中的模型参数进行求导，采用反向时间步进的隐式方法，逐步计算电场沿所有单元棱边相对模型参数的灵敏度转置，再结合系数插值矩阵，计算观测瞬变电磁响应数据对模型参数的灵敏度矩阵与向量的乘积；
模型更新参数的获取模块，用于结合观测瞬变电磁响应数据对模型参数的灵敏度，采用具有拟二次收敛速度的高速-牛顿法方程将目标函数转化为最小二乘问题，获取模型更新方向，并采用线性搜索获得最优模型更新步长；
模型更新模块，用于结合模型更新方向与最优模型更新步长计算模型参数，更新反演模型；
预测瞬变电磁响应数据的获取模块，用于基于矢量有限元方法对反演模型进行瞬变电磁三维正演，获取预测瞬变电磁响应数据；
误差评判模块，用于采用归一化误差评价预测瞬变电磁响应数据与观测瞬变电磁响应数据的拟合情况；若数据归一化拟合差达到预设阈值，反演方法终止；若未达到预设阈值，且当前迭代次数小于设定次数，则驱动观测瞬变电磁响应数据对模型参数的灵敏度计算模块执行。
进一步优选地，理论模型瞬变电磁三维正演单元包括：
有限元分析方程计算器，用于对时域的准静态Maxwell方程组两边取旋度，并采用基于非结构化四面体网格的一阶矢量有限元方法，获取有限元分析方程；其中，采用网格单元的棱边离散瞬变电磁发射源，使导线或线圈与网格单元的棱边重合；
有限元分析方程变形器，用于采用后退二阶欧拉方法对有限元分析方程中的时间进行离散，且通过泰勒级数展开，获取各时刻非均匀步长的电场时间导数表达式，对有限元分析方程进行变形；
条件设置器，用于采用Dirichlet边界条件，假设电场在理论模型的建模区域的边界上消失，设置电场为零的初始条件，且将理论模型的建模区域扩展，在测点和发射源的预设区域加密网格，其他区域逐步增大网格的尺寸，完成边界条件和初始条件的设定以及网格的设置；
电场计算器，用于基于边界条件和初始条件的设定以及网格的设置，采用直接求解器求解变形后有限元分析方程，计算各时刻的电场；
观测瞬变电磁响应数据计算器，用于根据电场和稀疏插值矩阵计算观测瞬变电磁响应数据。
进一步优选地，所述正则化目标函数为：
φ（m）=φd（m）+λφm（m）
φd（m）=0.5||Wd（F（m）-dobs）||2 L2
φm（m）=0.5||Wm（m-mref）||2 L2
其中，φd是数据拟合项；φm是模型稳定项；λ是正则化参数；m为反演模型参数；m=log（σ）；σ为电导率；Wd=diag（1/（dobs+ε））；F（m）表示初始反演模型参数m的正演响应，dobs为测点上的观测瞬变电磁响应数据，Wm为模型粗糙度矩阵，mref为含先验信息的参考模型参数；x是随机向量；J为观测瞬变电磁响应数据对模型参数的灵敏度；L2为二范数；ε为预设的较小值。
进一步优选地，所述电场沿所有单元棱边相对于模型参数的灵敏度为：
S n =∂En/∂m；
其中，S n 为第n个时间步电场沿所有单元棱边相对于模型参数的灵敏度；S n 的大小是N ed ×N m ，N ed 是网格棱边的数量；N m 是模型参数的数量；
Jn=QSn；
其中，Q是大小为N d×N ed的稀疏插值矩阵；Jn为第n个时间步观测瞬变电磁响应数据对模型参数的灵敏度；N d表示数据个数；数据个数为测取测点数据所用仪器对应时间道的数目与测点数目的乘积。
进一步优选地，第k次反演模型的模型参数为：
m k =m k-1+lδm k
其中，m k 为第k次反演模型的模型参数；m k-1为第k-1次反演模型的模型参数；δmk为模型更新方向；l为最优模型更新步长。
总体而言，通过本发明所构思的以上技术方案与现有技术相比，具有以下有益效果：
本发明提供了一种地面瞬变电磁法并行反演方法，在获取观测瞬变电磁响应数据时，对时域的准静态Maxwell方程组两边取旋度，采用基于非结构化四面体网格的一阶矢量有限元方法，获取有限元分析方程；采用Dirichlet边界条件，假设电场在建模区域的边界上消失，设置电场为零的初始条件，且将理论模型的建模区域扩展，在测点和发射源的预设区域加密网格，其他区域逐步增大网格的尺寸，完成边界条件和初始条件的设定以及网格的设置；可以精确模拟复杂地质结构的瞬变电磁响应。
本发明在正演过程中采用自适应时间步长技术，采用后退二阶欧拉方法对有限元分析方程中的时间进行离散，采用自适应时间步长方法，当以时间步进△t n-1计算一定次数的步长后，尝试把步长增加到△t n-1的k n 倍，且通过泰勒级数展开，获取各时刻非均匀步长的电场时间导数表达式，对有限元分析方程进行变形；并考虑到迭代法受限于方程的条件数，可能会出现不收敛的情况；本发明确定边界条件和初始条件后，可通过迭代法或直接法求解变形后的有限元分析方程得到某时刻的电场E n ，在求解过程中，当时间步长相同时，刚度矩阵保持不变，A n 有限元系数矩阵分解的结果可以被重复使用；由于每个△t n 的A n 的矩阵分解已经被存储因此，只需要执行直接求解器的回代求解阶段；在反演过程中，本发明使用反向时间步进的隐式方法求解灵敏度矩阵与向量之间的矩阵向量乘法时，矩阵分解结果A n 在反演过程中重新使用，结合MPI和OpenMP并行计算；综上所述，本发明极大地提高了正演与反演方法的效率。
本发明结合观测瞬变电磁响应数据对模型参数的灵敏度，采用具有拟二次收敛速度的高速-牛顿法方程将目标函数转化为最小二乘问题，具有拟二次收敛速度，获取模型更新方向，并采用线性搜索获得最优模型更新步长，具有良好的稳定性。
本发明利用伴随正演的方法计算观测瞬变电磁响应数据对模型参数的灵敏度矩阵与向量的乘积，极大降低了直接计算灵敏度矩阵的所需的巨大计算量和内存需求。因此，本发明所形成的反演方案可以高效率的解决复杂地质情况下的瞬变电磁三维反演问题。
附图说明
图1是本发明实施例提供的地面瞬变电磁法并行反演方法流程图.；
图2（a）是本发明实施例提供的全局网格对应的均匀半空间模型图，图中白色的框为核心区域；
图2（b）是本发明实施例提供的核心区域对应的均匀半空间模型图，发射框为发射源，大小为300m×300m，白色小点为测点，采用30m的观测点距，共81个测点；
图3（a）是本发明实施例提供的采用均匀半空间模型获取的坐标为（0，0）位置的测深曲线；
图3（b）是本发明实施例提供的图3（a）对应三维正演算法与解析解的相对误差示意图；
图3（c）是本发明实施例提供的采用均匀半空间模型获取的坐标为（-60m，-60m）位置的测深曲线；
图3（d）是本发明实施例提供的图3（c）对应三维正演算法与解析解的相对误差示意图；
图4（a）是本发明实施例提供的地面瞬变电磁探测装置，其中，采用大定回线源装置观测，发射回线框不动，发射电流等所有参数不变，接收部分（主机和接收线圈）在发射框内沿测线移动，测线与测点距离均为20米；
图4（b）是本发明实施例提供的地面瞬变电磁探测装置的发射波形，电流关断时间为2.58×10-4s；
图5是本发明实施例提供的测点分布情况，各方阵代表不同的发射源对应的测点，共有16个发射源；
图6（a）是本发明实施例提供的采用理论模型获取的三维视图；
图6（b）是本发明实施例提供的采用初始反演模型获取的三维视图；
图6（c）是本发明实施例提供的采用反演结果获取的穿过异常体中心切片的三维视图，其中，Y=700m；其中，白色虚线表示反演模型中两个不同层之间的界面，白色实线是两个异常体的轮廓；
图7是本发明实施例提供的反演收敛曲线；
图8是本发明实施例提供的反演结果预测数据的归一化拟合差分布图；
图9（a）~（c）是本发明实施例提供的随机选取三个测点的瞬变电磁测深曲线的数据拟合示意图。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。如图1所示，一方面，本发明提供了一种地面瞬变电磁法并行反演方法，包括以下步骤：
S1：获取用于反演的观测瞬变电磁响应数据；采集实测的实际观测瞬变电磁响应数据以及采用理论模型通过瞬变电磁三维正演获得理论观测瞬变电磁响应数据；需指出，这里的理论模型为现有的可以获取理论观测瞬变电磁响应数据的所有模拟模型；
其中，本发明中瞬变电磁三维正演采用非结构化四面体网格的矢量有限元方法；具体为：
从时域的准静态Maxwell方程组出发，忽略位移电流项得到：
▽×H=Je+Js；（1）
▽×E=-μ ∂H/∂t；（2）
其中，H为磁场强度；E为电场强度；μ为磁化率；Je为传导电流密度；Js为外加激励场源电流密度；
为了减小计算量，对Maxwell方程组两边取旋度，消去磁场项，可以得到如下扩散方程：
▽×▽×E+μσ ∂ Je（t）/∂t=-μ∂ Js（t）/∂t （3）
其中， σ为电导率；
为了求解方程（3），本发明采用基于非结构化四面体网格的一阶矢量有限元方法，获取有限元分析方程：
KE（t）+μMσ∂E（t）/∂t=-μM∂Js（t）/∂t（4）
其中，K为双旋度算子，Mσ和M为刚度矩阵，分别定义为：
其中，Ni和Nj是每个计算单元的域Ωe的基函数；Kij为双旋度算子的离散形式；Mij和 
为相关的质量矩阵；dv为计算单元体积；Ωe为计算单元域；方程（5）和（6）中的内积采用解析形式计算，方程（7）采用11个采样点的高斯求积计算；在公式（4）中，加载右边的源项（-μM∂ Js（t）/∂t）是至关重要的，在本发明中，使用网格单元的棱边离散瞬变电磁发射源，使导线或线圈与网格单元的棱边重合，发射源被分解成一系列电偶极子；方程（4）右边的源项只包含瞬变电磁发射源棱边上的值，其值与棱边长度和波形有关；
本发明采用后退二阶欧拉方法对有限元分析方程中的时间进行离散，为了减少所需模拟的时间步进次数，采用自适应时间步长方法；当以时间步进△t n-1计算一定次数的步长后，尝试把步长增加到△t n-1的k n 倍：
△t n =k n △t n-1 （8）
通过泰勒级数展开，在t n 时刻非均匀步长的电场时间导数可以表示为：
E n /∂t≈（（2kn+1）E n -（kn+1）2 E n-1+ kn 2 E n-2） /（（kn+1）△t n ）（9）
将方程（9）代入方程（4）得到：
公式（10）可以整理成如下形式：
A n E n =b n （11）
其中，A n 有限元系数矩阵为
；b n 为
；
求解方程（11）需要给定适当边界条件以指定建模域边界上的电磁场，以得到唯一的解，本发明采用Dirichlet边界条件，假设电场在理论模型的建模区域的边界上消失：
E（t）Ω=0 （12）
其中，Ω为建模域边界；
将建模区域扩展至一个数千米到数百千米的边界，以匹配Dirichlet边界条件的使用；边界的范围大小取决于发射源电极距的大小、所需模拟的最大时刻以及模型介质分布；但在测点和发射源的附近加密网格，以准确模拟地形和复杂三维地质体和确保计算精度；在核心区域以外逐步增大网格的尺寸，以减少网格数量提高效率；通过与均匀空间解析解对比确定边界范围、网格剖面情况是否合适；
求解方程（11）还需要合适的初始条件，本发明采用电场为零的初始条件，通过将离散的发射电流密度Js代入公式（11）实现发射源的加载，可以模拟任意发射波形；
确定边界条件和初始条件后，可通过迭代法或直接法求解方程（11）得到某时刻的电场E n ；考虑到迭代法受限于方程的条件数，可能会出现不收敛的情况；本发明中采用集群并行版本的Intel MKL PARDISO直接求解器求解方程（11），当时间步长相同时，刚度矩阵保持不变，A n 有限元系数矩阵分解的结果可以被重复使用；由于每个△t n 的A n 的矩阵分解已经被存储，因此只需要执行直接求解器的回代求解阶段；
此外，灵敏度矩阵（或其转置）与向量之间的矩阵-向量乘法涉及到向前和向后的时间步进；在这些过程中仍然需要用直接求解器求解类似方程组；
一旦通过求解方程（11）得到电场，测点的观测瞬变电磁数据（∂H/∂t）可以通过以下插值矩阵获得：
dobs=QE（13）
其中，Q是大小为N d×N ed的稀疏插值矩阵，N d和N ed表示数据个数和网格棱边数量；dobs为∂H/∂t；E为电场；数据个数为测点测取测点数据所用仪器对应时间道的数目与测点数目的乘积；
S2：利用非结构化四面体网格对反演空间进行剖分；其中，反演空间为测点及其附近预设延伸范围所对应的地下空间；
本发明采用非结构化四面体网格方法离散反演空间，通过使用非结构化网格，可以很容易地实现关键区域（如发射源和测点）的局部细化；发射源附近进行局部细化使得全场方法的实现更加容易，而测点周围的细化可以得到很精确的模拟结果；
S3：在S2中获取的非结构化四面体网格中设置电导率值，构建初始反演模型；
在三维反演中，初始反演模型的选择会影响反演效果和收敛速度，选择与真实模型相差较大的初始反演模型可能导致瞬变电磁反演不收敛，反演结果较差；因此，构建合适的初始反演模型对于获得高质量的反演结果具有重要意义；一般情况下可以参考已有的钻井、地质和地球物理先验信息建立简单的初始反演模型，但在先验信息较少、地质情况较复杂的地区，可以先进行一维反演，将一维反演结果作为初始模型；
S4：构建瞬变电磁正则化反演目标函数；
瞬变电磁数据三维反演是不适定问题，采用正则化目标函数可以获得具有地球物理意义的解；正则化目标函数为：
φ（m）=φd（m）+λφm（m）（14）
其中，φd是数据拟合项；φm是模型稳定项；λ是正则化参数；φd、φm和λ用于平衡数据拟合和模型约束；
φd和φm分别定义为：
φd（m）=0.5||Wd（F（m）-dobs）||2 L2 （15）
φm（m）=0.5||Wm（m-mref）||2 L2 （16）
其中，m为反演模型参数，为了保证反演过程中m为非负数，使得模型参数具有物理意义，取对数域的电导率作为模型参数，m=log（σ），Wd为数据加权矩阵，由于瞬变电磁数据跨越多个数量级，为了平衡不同时间通道的数据权重，采用测点观测瞬变电磁响应数据的导数进行构建，Wd=diag（1/（dobs+ε））；因为瞬变电磁数据可能会出现“变号”，因此在分母加上一个较小的值ε，可以避免数据权重出现无穷大值，ε可取最后一个时间道观测数据的平均数；F（m）表示初始反演模型参数m的正演响应，dobs为测点上的观测瞬变电磁响应数据，Wm为模型粗糙度矩阵，mref为含先验信息的参考模型参数；
正则化参数在权衡数据拟合项和模型约束项上起着重要作用，正则化参数的选择影响算法的稳定性和收敛速度；采用近似谱分析方法计算正则化参数：
其中，x是随机向量；在反演中，在初次迭代中使用公式（17）计算正则化参数；在后续迭代中，当收敛缓慢时使用冷却方法减小其值，或者重新使用公式（17）计算正则化参数；
根据距离反演单元最近的N个反演单元的体积与距离的权重构造粗糙度矩阵，表示如下：
其中，V i 是第i个反演单元的体积；V ij 表示第i个反演单元对应的第j个最近反演单元的体积；V ik 表示第i个反演单元对应的第k个最近反演单元的体积；r ij 表示第i个反演单元中心到第j个最近单元中心距离，通过以下公式计算：
其中，（x i ，y i ，z i ）和（x j ，y j ，z j ）为第i个反演和第j个反演单元中心的中心坐标；当在上述方法中，可以通过参数N控制粗糙度矩阵反演结果的平滑度，但N越大时反演结果越光滑；
S5：利用伴随正演的方法计算灵敏度矩阵（或灵敏度矩阵的转置）与向量的乘积；
灵敏度矩阵表示电场沿所有单元棱边相对模型的变化量的变化，对方程（10）和（11）中模型参数m求导可以得到：
定义电场沿所有单元棱边相对于模型参数的灵敏度为S，S n 表示第n个时间步的灵敏度：
S n =∂En/∂m（21）
其中，S n 的大小是N ed ×N m ，N ed 是棱边的数量；N m 是模型参数的数量；
为了简洁表示，还作如下定义：
将公式（21）、公式（22）和公式（23）代入公式（20）得到：
其中：S n-1和S n-2分别表示第n-1和n-2个时间步的电场灵敏度；矩阵G n的列数取决于模型参数的个数N m ，N m 通常很大（现实问题从数万到数百万）；因此，直接求解方程（24）的计算非常大；
将观测瞬变电磁响应数据对模型参数的灵敏度定义为J，Jn第n个时间步的观测数据灵敏度，大小为N d ×N m ，N d 是观测数据的数量，N m 是模型参数的数量；Sn和Jn的关系通过下式表示：
Jn=QSn（25）
数据在所有时间道的灵敏度可以写成一个整体：
通常数据的数量远小于模型参数的数量；因此，取公式（26）的转置，得到灵敏度矩阵的转置：
其中，T代表转置；
对公式（27）整体进行矩阵因式分解成本极高；因此，通常采取逐步求解灵敏度矩阵的转置；与正演问题求解方程（11）的正向步进的方式相反，使用反向步进方法求解灵敏度；公式（27）右端项依赖于列数更少的QT；然而，对于实际瞬变电磁反演问题，灵敏度矩阵或其转置的显式计算仍然是非常困难；本发明使用反向时间步进的隐式方法求解灵敏度矩阵时，采用计算灵敏度矩阵（或其转置）与向量之间的矩阵向量乘法；在正演模拟中，采用自适应时间步进方法只需对矩阵An进行几次更新，其矩阵分解被存储起来；在反演过程中，当计算J（或JT）和向量之间的矩阵向量乘法时，这些矩阵分解结果将在反演过程中重新使用，每次计算J矩阵与向量的乘积相当于一次正演计算，结合MPI和OpenMP并行计算可以有效加速反演；
S6：采用LSQR方法求解等价高斯-牛顿法方程的最小二乘问题，获得模型更新方向，线性搜索获得最优模型更新步长，对预测模型进行更新；
采用具有拟二次收敛速度的高斯-牛顿方法对目标函数进行最小化，将泰勒展开公式应用于目标函数，并忽略高阶项可以获得以下法方程：
其中，δmk为模型更新方向，数据误差r k-1 =F（mk-1）-dobs，k为迭代次数；将法方程转换为如下最小二乘形式：
采用最小二乘（LSQR）算法求解方程（29）获得模型更新量，采用线性搜索方法找到优化的模型更新步长l，模型更新可以表示为：
m k =m k-1+lδm k （30）
S7：基于矢量有限元法对预测模型m k 进行瞬变电磁三维正演，得到预测模型瞬变电磁响应数据；在正演中结合了自适应时间步长技术、Intel MKL Pardiso直接求解器、MPI和OpenMP并行计算技术提高计算效率；
S8：采用归一化误差评价模型预测数据与观测数据的拟合情况，判断是否达到预设阈值，判断是否达到最大迭代次数；
判断观测瞬变电磁响应数据和预测瞬变电磁响应数据是否拟合，观测瞬变电磁响应数据与预测瞬变电磁响应数据d pre之间的数据归一化拟合差可用以下公式表示：
数据归一化拟合差达到设定阈值或收敛停止时，反演算法终止，输出反演结果；若未达到设定的最小数据归一化拟合差，则判断是否达到最大的反演迭代次数，若当前迭代次数小于设定次数，则循环进行下一次迭代，若到达设定次数，则停止反演迭代，输出反演结果；
S9：迭代执行S5~S8直至到达设定的终止条件推出，获得电导率模型，实现瞬变电磁三维反演。
另一方面，本发明提供了一种地面瞬变电磁法并行反演系统，其特征在于，包括：
瞬变电磁数据的获取模块，用于获得用于反演的观测瞬变电磁响应数据，包括野外实地数据采集单元与理论模型瞬变电磁三维正演单元，分别获取实际观测瞬变电磁响应数据和理论观测瞬变电磁响应数据；
空间剖分模块，用于将反演空间剖分成非结构化四面体网格；其中，反演空间为测点及其附近预设延伸范围所对应的地下空间；
初始反演模型的构建模块，用于在非结构化四面体网格中设置电导率值，构建初始反演模型，并结合观测瞬变电磁响应数据构建初始反演模型的正则化目标函数；
观测瞬变电磁响应数据对模型参数的灵敏度计算模块，用于对棱边电场方程中初始反演模型中的模型参数进行求导，采用反向时间步进的隐式方法，逐步计算电场沿所有单元棱边相对模型参数的灵敏度转置，再结合系数插值矩阵，计算观测瞬变电磁响应数据对模型参数的灵敏度矩阵与向量的乘积；
模型更新参数的获取模块，用于结合观测瞬变电磁响应数据对模型参数的灵敏度，采用具有拟二次收敛速度的高速-牛顿法方程将目标函数转化为最小二乘问题，获取模型更新方向，并采用线性搜索获得最优模型更新步长；
模型更新模块，用于结合模型更新方向与最优模型更新步长计算模型参数，更新反演模型；
预测瞬变电磁响应数据的获取模块，用于基于矢量有限元方法对反演模型进行瞬变电磁三维正演，获取预测瞬变电磁响应数据；
误差评判模块，用于采用归一化误差评价预测瞬变电磁响应数据与观测瞬变电磁响应数据的拟合情况；若数据归一化拟合差达到预设阈值，反演方法终止；若未达到预设阈值，且当前迭代次数小于设定次数，则驱动观测瞬变电磁响应数据对模型参数的灵敏度计算模块执行。
进一步优选地，所述瞬变电磁数据的获取模块包括野外实地数据采集单元和观测瞬变电磁响应数据获取单元；
其中，理论模型瞬变电磁三维正演单元包括：
有限元分析方程计算器，用于对时域的准静态Maxwell方程组两边取旋度，并采用基于非结构化四面体网格的一阶矢量有限元方法，获取有限元分析方程；其中，采用网格单元的棱边离散瞬变电磁发射源，使导线或线圈与网格单元的棱边重合；
有限元分析方程变形器，用于采用后退二阶欧拉方法对有限元分析方程中的时间进行离散，且通过泰勒级数展开，获取各时刻非均匀步长的电场时间导数表达式，对有限元分析方程进行变形；
条件设置器，用于采用Dirichlet边界条件，假设电场在建模区域的边界上消失，设置电场为零的初始条件，且将理论模型的建模区域扩展，在测点和发射源的预设区域加密网格，其他区域逐步增大网格的尺寸，完成边界条件和初始条件的设定以及网格的设置；
电场计算器，用于基于边界条件和初始条件的设定以及网格的设置，采用直接求解器求解变形后有限元分析方程，计算各时刻的电场；
观测瞬变电磁响应数据计算器，用于根据电场和稀疏插值矩阵计算观测瞬变电磁响应数据。
进一步优选地，所述正则化目标函数为：
φ（m）=φd（m）+λφm（m）
φd（m）=0.5||Wd（F（m）-dobs）||2 L2
φm（m）=0.5||Wm（m-mref）||2 L2
其中，φd是数据拟合项；φm是模型稳定项；λ是正则化参数；m为反演模型参数；m=log（σ）；σ为电导率；Wd=diag（1/（dobs+ε））；F（m）表示初始反演模型参数m的正演响应，dobs为测点上的观测瞬变电磁响应数据，Wm为模型粗糙度矩阵，mref为含先验信息的参考模型参数；x是随机向量；J为灵敏度矩阵；L2为二范数；ε为预设的较小值。
进一步优选地，所述电场沿所有单元棱边相对于模型参数的灵敏度为：
S n =∂En/∂m；
其中，S n 为第n个时间步电场沿所有单元棱边相对于模型参数的灵敏度；S n 的大小是N ed ×N m ，N ed 是网格棱边的数量；N m 是模型参数的数量；
Jn=QSn；
其中，Q是大小为N d×N ed的稀疏插值矩阵；J为观测瞬变电磁响应数据对模型参数的灵敏度；N d表示数据个数；数据个数为测取测点数据所用仪器对应时间道的数目与测点数目的乘积。
进一步优选地，第k次反演模型的模型参数为：
m k =m k-1+lδm k
其中，m k 为第k次反演模型的模型参数；m k-1为第k-1次反演模型的模型参数；δm k 为模型更新方向；l为最优模型更新步长。
实施例1
为了验证本发明所提供的瞬变电磁三维正演的有效性，本实施例提供了一个均匀半空间模型如图2（a）和图2（b）所示，地下空间电导率为 0.1S/m，空气层电导率为
S/m；图 2（b）中的白色线框为发射线框为发射源，大小为 300×300m，白色为测点，采用 30m 的观测点距，共 81 个测点；发射波形为斜阶跃波，电流关断时间为 5×10-6 s，观测时间取等对数间隔 10-4~10-1.5 s共 25 个时间道；数据数目为 81×25=2025 个；
本实施例采用解析法获取解析解，在高性能计算集群上完成，每个节点包含2个Inter(R) Xeon(R) 2.5GHz的 CPU，每个 CPU 包含 20 个核，单节点内存为 384GB；使用1个节点进行计算，采用自适应时间步进方法正演时，斜阶跃波形被离散为 10 段，初始时间步长为5×10-7 s，每计算20次等步长 
后尝试将时间步长变为当前步长 2 倍，通过241 个步进后完成计算，共有12个不同步长，计算时间为154 s；采用等步长计算时，需要计算 63246个步进，约需要4小时50分钟；
该理论模型使用本发明所提供的正演方法得到的理论观测瞬变电磁响应数据与解析法得到的解析解对比；如图3（a）~（d）所示，瞬变电磁三维正演与解析解的相对误差非常小，大部分时间道上数据相对误差小于2%；研究结果验证了自适应步长方法在确保精度的基础上极大提高了正演计算效率，可以满足实际瞬变电磁数据解释的需求，为瞬变电磁三维反演提供了良好的基础。
实施例2
根据图4（a）与图5地面瞬变电测观测装置建立了如图6（a）所示的理论模型；该理论模型主要由上下两层构成，从而模拟该测区覆盖层与基岩总体地质情况，电导率分别为0.01S/m和0.002S/m。第一层设置了一个高阻异常体和一个低阻的水平板状，异常体的电阻率分布为0.001S/m和0.2S/m；该理论模型采用如图4（b）所示的发射波形，电流关断时间为2.58×10-4 s；
根据图4（a）与图5地面瞬变电测观测装置建立了如图6（b）所示的反演空间，并利用非结构化四面体网格离散反演空间；反演区域X方向为-500m~2500m，Y方向为-240m~1660m，从地形表面到Z = 250m(地形的最高高程为1468m)。整个模型区域大小为10km×10km×0km，X方向从-4000~6000m，Y方向-4500m~5500m，Z方向为-4000m~6000m。整个模型的网格包含715636个单元和837234条棱边，反演域有417020个单元；如图6（b）所示，图6（a）理论模型的反演中初始模型为电导率为0.01S/m的均匀空间；
通过S4~S9，图6（a）理论模型结果25次反演迭代后，反演结果的三维视图如图6（c）所示；从图中可以看出，三维反演对两个异常体和双层结构得到恢复；理论模型的反演结果验证了本发明能够对复杂地形地区的瞬变电磁数据进行精细解释；
反演收敛曲线如图7所示，从图中可以看到，经过25次高斯牛顿迭代后，归一化拟合收敛到2%以下；在一个集群节点中使用2个MPI进程和20个OpenMP线程，计算时间约为67个小时；对于该模型每一次高斯-牛顿迭代中使用10次LSQR迭代来进行模型更新；
反演结果预测数据的归一化拟合差如图8所示，所有测点的归一化拟合差小于5%；
部分测点的瞬变电磁测深曲线的数据拟合如图9（a）~9（c）所示，从图中可以看到测点的瞬变电磁测深曲线的数据拟合较好。
综上所述，本发明与现有技术相比，存在以下优势：
本发明提供了一种地面瞬变电磁法并行反演方法，在获取观测瞬变电磁响应数据时，对时域的准静态Maxwell方程组两边取旋度，采用基于非结构化四面体网格的一阶矢量有限元方法，获取有限元分析方程；采用Dirichlet边界条件，假设电场在建模区域的边界上消失，设置电场为零的初始条件，且将理论模型的建模区域扩展，在测点和发射源的预设区域加密网格，其他区域逐步增大网格的尺寸，完成边界条件和初始条件的设定以及网格的设置；可以精确模拟复杂地质结构的瞬变电磁响应。
本发明在正演过程中采用自适应时间步长技术，采用后退二阶欧拉方法对有限元分析方程中的时间进行离散，采用自适应时间步长方法，当以时间步进△t n-1计算一定次数的步长后，尝试把步长增加到△t n-1的k n 倍，且通过泰勒级数展开，获取各时刻非均匀步长的电场时间导数表达式，对有限元分析方程进行变形；并考虑到迭代法受限于方程的条件数，可能会出现不收敛的情况；本发明确定边界条件和初始条件后，可通过迭代法或直接法求解变形后的有限元分析方程得到某时刻的电场E n ，在求解过程中，当时间步长相同时，刚度矩阵保持不变，A n 有限元系数矩阵分解的结果可以被重复使用；由于每个△t n 的A n 的矩阵分解已经被存储，因此只需要执行直接求解器的回代求解阶段；在反演过程中，本发明使用反向时间步进的隐式方法求解灵敏度矩阵与向量之间的矩阵向量乘法时，矩阵分解结果A n 在反演过程中重新使用，结合MPI和OpenMP并行计算；综上所述，本发明极大地提高了正演与反演方法的效率。
本发明结合观测瞬变电磁响应数据对模型参数的灵敏度，采用具有拟二次收敛速度的高速-牛顿法方程将目标函数转化为最小二乘问题，具有拟二次收敛速度，获取模型更新方向，并采用线性搜索获得最优模型更新步长，具有良好的稳定性。
本发明利用伴随正演的方法计算观测瞬变电磁响应数据对模型参数的灵敏度矩阵与向量的乘积，极大降低了直接计算灵敏度矩阵的所需的巨大计算量和内存需求。因此，本发明所形成的反演方案可以高效率的解决复杂地质情况下的瞬变电磁三维反演问题。
本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。
Claims (10)
1.一种地面瞬变电磁法并行反演方法，其特征在于，包括以下步骤：
S1：获取用于反演的观测瞬变电磁响应数据，包括采集地面的实际观测瞬变电磁响应数据和采用理论模型通过瞬变电磁三维正演获取理论观测瞬变电磁响应数据；其中，理论模型为能够获取理论观测瞬变电磁响应数据的模拟模型；
S2：将反演空间剖分成非结构化四面体网格；其中，反演空间为测点及其附近预设延伸范围所对应的地下空间；
S3：在非结构化四面体网格中设置电导率值，构建初始反演模型，并结合观测瞬变电磁响应数据构建初始反演模型的正则化目标函数；
S4：对棱边电场方程中初始反演模型中的模型参数进行求导，采用反向时间步进的隐式方法，逐步计算电场沿所有单元棱边相对模型参数的灵敏度转置，再结合系数插值矩阵，计算观测瞬变电磁响应数据对模型参数的灵敏度矩阵与向量的乘积；
S5：结合观测瞬变电磁响应数据对模型参数的灵敏度与向量的乘积，采用具有拟二次收敛速度的高速-牛顿法方程将目标函数转化为最小二乘问题，获取模型更新方向，并采用线性搜索获得最优模型更新步长；
S6：结合模型更新方向与最优模型更新步长计算模型参数，更新反演模型；
S7：基于矢量有限元方法对更新后的反演模型进行瞬变电磁三维正演，获取预测瞬变电磁响应数据；
S8：采用归一化误差评价预测瞬变电磁响应数据与观测瞬变电磁响应数据的拟合情况；若数据归一化拟合差达到预设阈值，反演方法终止；若未达到预设阈值，且当前迭代次数小于设定次数，则转至S4。
2.根据权利要求1所述的地面瞬变电磁法并行反演方法，其特征在于，获取所述理论观测瞬变电磁响应数据的方法，包括以下步骤：
对时域的准静态Maxwell方程组两边取旋度，并采用基于非结构化四面体网格的一阶矢量有限元方法，获取有限元分析方程；其中，采用网格单元的棱边离散瞬变电磁发射源，使导线或线圈与网格单元的棱边重合；
采用后退二阶欧拉方法对有限元分析方程中的时间进行离散，且通过泰勒级数展开，获取各时刻非均匀步长的电场时间导数表达式，对有限元分析方程进行变形；
采用Dirichlet边界条件，假设电场在理论模型的建模区域边界上消失，设置电场为零的初始条件，且将理论模型的建模区域扩展，在测点和发射源的预设区域加密网格，其他区域逐步增大网格的尺寸，完成边界条件和初始条件的设定以及网格的设置；
基于边界条件和初始条件的设定以及网格的设置，采用直接求解器求解变形后有限元分析方程，计算各时刻的电场；
根据电场和稀疏插值矩阵计算理论观测瞬变电磁响应数据。
3.根据权利要求1或2所述的地面瞬变电磁法并行反演方法，其特征在于，所述正则化目标函数为：
φ（m）=φd（m）+λφm（m）
φd（m）=0.5||Wd（F（m）-dobs）||2 L2
φm（m）=0.5||Wm（m-mref）||2 L2
其中，φd是数据拟合项；φm是模型稳定项；λ是正则化参数；m为反演模型参数；m=log（σ）；σ为电导率；Wd=diag（1/（dobs+ε））；F（m）表示初始反演模型参数m的正演响应，dobs为测点上的观测瞬变电磁响应数据，Wm为模型粗糙度矩阵，mref为含先验信息的参考模型参数；x是随机向量；J为观测瞬变电磁响应数据对模型参数的灵敏度；L2为二范数；ε为预设的较小值。
4.根据权利要求3所述的地面瞬变电磁法并行反演方法，其特征在于，所述电场沿所有单元棱边相对于模型参数的灵敏度为：
Sn=∂En/∂m；
其中，Sn为第n个时间步电场沿所有单元棱边相对于模型参数的灵敏度；S n的大小是N ed ×N m ，N ed 是网格棱边的数量；N m 是模型参数的数量；
Jn=QSn；
其中，Q是大小为N d×N ed的稀疏插值矩阵；Jn为第n个时间步观测瞬变电磁响应数据对模型参数的灵敏度；N d表示数据个数；数据个数为测取测点数据所用仪器对应时间道的数目与测点数目的乘积。
5.根据权利要求4所述的地面瞬变电磁法并行反演方法，其特征在于，第k次反演模型的模型参数为：
m k =m k-1+lδm k
其中，m k 为第k次反演模型的模型参数；m k-1为第k-1次反演模型的模型参数；δm k 为模型更新方向；l为最优模型更新步长。
6.一种地面瞬变电磁法并行反演系统，其特征在于，包括：
瞬变电磁数据的获取模块，用于获得用于反演的观测瞬变电磁响应数据，包括野外实地数据采集单元与理论模型瞬变电磁三维正演单元，分别获取实际观测瞬变电磁响应数据和理论观测瞬变电磁响应数据；理论模型为能够获取理论观测瞬变电磁响应数据的模拟模型；
空间剖分模块，用于将反演空间剖分成非结构化四面体网格；其中，反演空间为测点及其附近预设延伸范围所对应的地下空间；
初始反演模型的构建模块，用于在非结构化四面体网格中设置电导率值，构建初始反演模型，并结合观测瞬变电磁响应数据构建初始反演模型的正则化目标函数；
观测瞬变电磁响应数据对模型参数的灵敏度计算模块，用于对棱边电场方程中初始反演模型中的模型参数进行求导，采用反向时间步进的隐式方法，逐步计算电场沿所有单元棱边相对模型参数的灵敏度转置，再结合系数插值矩阵，计算观测瞬变电磁响应数据对模型参数的灵敏度矩阵与向量的乘积；
模型更新参数的获取模块，用于结合观测瞬变电磁响应数据对模型参数的灵敏度与向量的乘积，采用具有拟二次收敛速度的高速-牛顿法方程将目标函数转化为最小二乘问题，获取模型更新方向，并采用线性搜索获得最优模型更新步长；
模型更新模块，用于结合模型更新方向与最优模型更新步长计算模型参数，更新反演模型；
预测瞬变电磁响应数据的获取模块，用于基于矢量有限元方法对更新后的反演模型进行瞬变电磁三维正演，获取预测瞬变电磁响应数据；
误差评判模块，用于采用归一化误差评价预测瞬变电磁响应数据与观测瞬变电磁响应数据的拟合情况；若数据归一化拟合差达到预设阈值，反演方法终止；若未达到预设阈值，且当前迭代次数小于设定次数，则驱动观测瞬变电磁响应数据对模型参数的灵敏度计算模块执行。
7.根据权利要求6所述的地面瞬变电磁法并行反演系统，其特征在于，所述理论模型瞬变电磁三维正演单元包括：
有限元分析方程计算器，用于对时域的准静态Maxwell方程组两边取旋度，并采用基于非结构化四面体网格的一阶矢量有限元方法，获取有限元分析方程；其中，采用网格单元的棱边离散瞬变电磁发射源，使导线或线圈与网格单元的棱边重合；
有限元分析方程变形器，用于采用后退二阶欧拉方法对有限元分析方程中的时间进行离散，且通过泰勒级数展开，获取各时刻非均匀步长的电场时间导数表达式，对有限元分析方程进行变形；
条件设置器，用于采用Dirichlet边界条件，假设电场在理论模型的建模区域边界上消失，设置电场为零的初始条件，且将理论模型的建模区域扩展，在测点和发射源的预设区域加密网格，其他区域逐步增大网格的尺寸，完成边界条件和初始条件的设定以及网格的设置；
电场计算器，用于基于边界条件和初始条件的设定以及网格的设置，采用直接求解器求解变形后有限元分析方程，计算各时刻的电场；
观测瞬变电磁响应数据计算器，用于根据电场和稀疏插值矩阵计算理论观测瞬变电磁响应数据。
8.根据权利要求7所述的地面瞬变电磁法并行反演系统，其特征在于，所述正则化目标函数为：
φ（m）=φd（m）+λφm（m）
φd（m）=0.5||Wd（F（m）-dobs）||2 L2
φm（m）=0.5||Wm（m-mref）||2 L2
其中，φd是数据拟合项；φm是模型稳定项；λ是正则化参数；m为反演模型参数；m=log（σ）；σ为电导率；Wd=diag（1/（dobs+ε））；F（m）表示初始反演模型参数m的正演响应， dobs为测点上的观测瞬变电磁响应数据，Wm为模型粗糙度矩阵，mref为含先验信息的参考模型参数；x是随机向量；J为观测瞬变电磁响应数据对模型参数的灵敏度；L2为二范数；ε为预设的较小值。
9.根据权利要求8所述的地面瞬变电磁法并行反演系统，其特征在于，所述电场沿所有单元棱边相对于模型参数的灵敏度为：
Sn=∂En/∂m；
其中，Sn为第n个时间步电场沿所有单元棱边相对于模型参数的灵敏度；S n的大小是N ed ×N m ，N ed 是网格棱边的数量；N m 是模型参数的数量；
Jn=QSn；
其中，Q是大小为N d×N ed的稀疏插值矩阵；Jn为第n个时间步观测瞬变电磁响应数据对模型参数的灵敏度；N d表示数据个数；数据个数为测取测点数据所用仪器对应时间道的数目与测点数目的乘积。
10.根据权利要求9所述的地面瞬变电磁法并行反演系统，其特征在于，第k次反演模型的模型参数为：
m k =m k-1+lδm k
其中，m k 为第k次反演模型的模型参数；m k-1为第k-1次反演模型的模型参数；δm k 为模型更新方向；l为最优模型更新步长。
Priority Applications (2)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210888762.4A CN114970290B (zh) | 2022-07-27 | 2022-07-27 | 一种地面瞬变电磁法并行反演方法及系统 |
| US18/088,630 US20240054265A1 (en) | 2022-07-27 | 2022-12-26 | Parallel inversion method and system for ground-based transient electromagnetic method |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210888762.4A CN114970290B (zh) | 2022-07-27 | 2022-07-27 | 一种地面瞬变电磁法并行反演方法及系统 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN114970290A true CN114970290A (zh) | 2022-08-30 |
| CN114970290B CN114970290B (zh) | 2022-11-01 |
Family
ID=82970091
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202210888762.4A Active CN114970290B (zh) | 2022-07-27 | 2022-07-27 | 一种地面瞬变电磁法并行反演方法及系统 |
Country Status (2)
| Country | Link |
|---|---|
| US (1) | US20240054265A1 (zh) |
| CN (1) | CN114970290B (zh) |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117436301A (zh) * | 2023-09-21 | 2024-01-23 | 长江大学 | 基于非结构有限元考虑激电效应的时移大地电磁正演方法 |
| CN117454675A (zh) * | 2023-12-26 | 2024-01-26 | 中国地质科学院地球物理地球化学勘查研究所 | 一种定源瞬变电磁响应校正方法、系统及设备 |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104537714A (zh) * | 2015-01-07 | 2015-04-22 | 吉林大学 | 磁共振与瞬变电磁空间约束联合反演方法 |
| US20190129057A1 (en) * | 2017-10-30 | 2019-05-02 | Baker Hughes, A Ge Company, Llc | Multiple casing inspection tool combination with 3d arrays and adaptive dual operational modes |
| CN112949134A (zh) * | 2021-03-09 | 2021-06-11 | 吉林大学 | 基于非结构有限元方法的地-井瞬变电磁反演方法 |
| US20210357540A1 (en) * | 2020-05-15 | 2021-11-18 | International Business Machines Corporation | Matrix sketching using analog crossbar architectures |
| CN114779355A (zh) * | 2022-02-24 | 2022-07-22 | 中国地质大学（武汉） | 基于发射电流全波形的地面瞬变电磁法反演方法及装置 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| KR101092668B1 (ko) * | 2009-06-17 | 2011-12-13 | 서울대학교산학협력단 | 파형 역산을 이용한 지하 구조의 영상화 장치와 방법 |
-
2022
- 2022-07-27 CN CN202210888762.4A patent/CN114970290B/zh active Active
- 2022-12-26 US US18/088,630 patent/US20240054265A1/en not_active Abandoned
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104537714A (zh) * | 2015-01-07 | 2015-04-22 | 吉林大学 | 磁共振与瞬变电磁空间约束联合反演方法 |
| US20190129057A1 (en) * | 2017-10-30 | 2019-05-02 | Baker Hughes, A Ge Company, Llc | Multiple casing inspection tool combination with 3d arrays and adaptive dual operational modes |
| US20210357540A1 (en) * | 2020-05-15 | 2021-11-18 | International Business Machines Corporation | Matrix sketching using analog crossbar architectures |
| CN112949134A (zh) * | 2021-03-09 | 2021-06-11 | 吉林大学 | 基于非结构有限元方法的地-井瞬变电磁反演方法 |
| CN114779355A (zh) * | 2022-02-24 | 2022-07-22 | 中国地质大学（武汉） | 基于发射电流全波形的地面瞬变电磁法反演方法及装置 |
Non-Patent Citations (2)
| Title |
|---|
| XIAOYUE CAO 等: "3-D Inversion of Z-Axis Tipper Electromagnetic Data Using Finite-Element Method With Unstructured Tetrahedral Grids", 《IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING ( VOLUME: 60)》 * |
| 胡祥云 等: "大地电磁三维数据空间反演并行算法研究", 《地球物理学报》 * |
Cited By (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117436301A (zh) * | 2023-09-21 | 2024-01-23 | 长江大学 | 基于非结构有限元考虑激电效应的时移大地电磁正演方法 |
| CN117436301B (zh) * | 2023-09-21 | 2024-04-16 | 长江大学 | 基于非结构有限元考虑激电效应的时移大地电磁正演方法 |
| CN117454675A (zh) * | 2023-12-26 | 2024-01-26 | 中国地质科学院地球物理地球化学勘查研究所 | 一种定源瞬变电磁响应校正方法、系统及设备 |
| CN117454675B (zh) * | 2023-12-26 | 2024-04-09 | 中国地质科学院地球物理地球化学勘查研究所 | 一种定源瞬变电磁响应校正方法、系统及设备 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN114970290B (zh) | 2022-11-01 |
| US20240054265A1 (en) | 2024-02-15 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN114970290B (zh) | 一种地面瞬变电磁法并行反演方法及系统 | |
| Hunt et al. | Efficient data assimilation for spatiotemporal chaos: A local ensemble transform Kalman filter | |
| CN114966685B (zh) | 基于InSAR和深度学习的大坝形变监测及预测方法 | |
| CN112949134B (zh) | 基于非结构有限元方法的地-井瞬变电磁反演方法 | |
| Liu et al. | Model updating of complex structures using the combination of component mode synthesis and Kriging predictor | |
| Gillman et al. | A fast direct solver for quasi-periodic scattering problems | |
| Le Bouteiller et al. | A discontinuous Galerkin fast-sweeping eikonal solver for fast and accurate traveltime computation in 3D tilted anisotropic media | |
| CN115238550A (zh) | 自适应非结构网格的滑坡降雨的地电场数值模拟计算方法 | |
| Xu et al. | An efficient method for statistical moments and reliability assessment of structures | |
| CN109490948A (zh) | 地震声学波动方程矢量并行计算方法 | |
| Lin et al. | Three-dimensional conjugate gradient inversion of magnetotelluric sounding data | |
| Luo | A POD-based reduced-order TSCFE extrapolation iterative format for two-dimensional heat equations | |
| Wang et al. | Geophysical electromagnetic modeling and evaluation: a review | |
| Likanapaisal et al. | Dynamic data integration and quantification of prediction uncertainty using statistical-moment equations | |
| Zhang et al. | A new 3-D ray tracing method based on LTI using successive partitioning of cell interfaces and traveltime gradients | |
| Karimaghaei et al. | Boundary integral formulation of the standard eigenvalue problem for the 2-D Helmholtz equation | |
| Li et al. | A frozen Gaussian approximation-based multi-level particle swarm optimization for seismic inversion | |
| CN118033764B (zh) | 一种多尺度多方法地质体电阻率三维成像方法及系统 | |
| CN110794469A (zh) | 基于最小地质特征单元约束的重力反演方法 | |
| Chan et al. | Reduced storage nodal discontinuous Galerkin methods on semi-structured prismatic meshes | |
| Nintcheu Fata et al. | A fast spectral Galerkin method for hypersingular boundary integral equations in potential theory | |
| BR102016024241B1 (pt) | Método para modelar um sistema geofísico e método para determinar propriedades físicas | |
| Bergamaschi et al. | Spectral acceleration of parallel iterative eigensolvers for large scale scientific computing | |
| CN118033764A (zh) | 一种多尺度多方法地质体电阻率三维成像方法及系统 | |
| Liu et al. | Three-Dimensional Inversion of Time-Domain Electromagnetic Data Using Various Loop Source Configurations |
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 | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |