CN115508898A - G-s变换的接地长导线源瞬变电磁快速正反演方法及系统 - Google Patents

G-s变换的接地长导线源瞬变电磁快速正反演方法及系统 Download PDF

Info

Publication number
CN115508898A
CN115508898A CN202211243244.3A CN202211243244A CN115508898A CN 115508898 A CN115508898 A CN 115508898A CN 202211243244 A CN202211243244 A CN 202211243244A CN 115508898 A CN115508898 A CN 115508898A
Authority
CN
China
Prior art keywords
formula
inversion
magnetic field
source
follows
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
Application number
CN202211243244.3A
Other languages
English (en)
Other versions
CN115508898B (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.)
Sichuan Zhongshui Chengkanyuan Geophysical Exploration Co ltd
Original Assignee
Sichuan Zhongshui Chengkanyuan Geophysical Exploration Co ltd
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 Sichuan Zhongshui Chengkanyuan Geophysical Exploration Co ltd filed Critical Sichuan Zhongshui Chengkanyuan Geophysical Exploration Co ltd
Priority to CN202211243244.3A priority Critical patent/CN115508898B/zh
Publication of CN115508898A publication Critical patent/CN115508898A/zh
Application granted granted Critical
Publication of CN115508898B publication Critical patent/CN115508898B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

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
    • 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
    • G01V3/10—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 using induction coils
    • 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
    • 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/14—Fourier, Walsh or analogous domain transformations, e.g. Laplace, Hilbert, Karhunen-Loeve, transforms
    • 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
    • G06F9/00—Arrangements for program control, e.g. control units
    • G06F9/06—Arrangements for program control, e.g. control units using stored programs, i.e. using an internal store of processing equipment to receive or retain programs
    • G06F9/46—Multiprogramming arrangements
    • G06F9/50—Allocation of resources, e.g. of the central processing unit [CPU]
    • G06F9/5005—Allocation of resources, e.g. of the central processing unit [CPU] to service a request
    • G06F9/5027—Allocation of resources, e.g. of the central processing unit [CPU] to service a request the resource being a machine, e.g. CPUs, Servers, Terminals
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2209/00—Indexing scheme relating to G06F9/00
    • G06F2209/50—Indexing scheme relating to G06F9/50
    • G06F2209/5018—Thread allocation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Data Mining & Analysis (AREA)
  • Pure & Applied Mathematics (AREA)
  • Computational Mathematics (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Databases & Information Systems (AREA)
  • Algebra (AREA)
  • Remote Sensing (AREA)
  • Geophysics (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • Geology (AREA)
  • Environmental & Geological Engineering (AREA)
  • Operations Research (AREA)
  • Electromagnetism (AREA)
  • Computing Systems (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本发明涉及G‑S变换的接地长导线源瞬变电磁快速正反演方法及系统，属于地球物理电磁勘探技术领域，针对地面接地长导线源瞬变电磁法，提供一种基于G‑S变换的接地长导线源瞬变电磁快速正反演方法技术流程，设计开发了一套基于G‑S变换的接地长导线源瞬变电磁法快速正反演解释软件系统。主要通过包含OMP并行计算的G‑S算法加快正演计算过程，通过单纯的OMP并行计算加速反演过程，极大的提高了接地长导线源瞬变电磁法正反演计算的效率，解决了反演时间长、效率低的问题。为快速获得现场接地长导线源瞬变电磁勘探数据反演结果提供了技术支撑。

Description

G-S变换的接地长导线源瞬变电磁快速正反演方法及系统
技术领域
本发明属于地球物理电磁勘探技术领域，具体涉及G-S变换的接地长导线源瞬变电磁快速正反演方法及系统。本发明主要解决了接地长导线源瞬变电磁法反演时间长、效率低的问题。
背景技术
接地长导线源瞬变电磁勘探方法是一种探测地下目标层位或者异常体的位置和赋存状态的地球物理勘探方法。主要利用接地导线源作为激发设备，向地下激发阶跃电流信号从而在空间产生一次磁场，一次磁场通过激励地下地质体感应出变化的二次场，然后在远离导线源的位置利用接收线圈采集二次场信号。这些二次场信号中包含的目标地质体信息，利用地球物理正反演方法可以提取这些地质信息，达到探测地下目标体的目的。通常情况下，接地长导线源瞬变电磁勘探方法由于常规正反演计算时间长，在野外现场工作时无法及时给出正反演结果。
因此，现阶段需设计G-S变换的接地长导线源瞬变电磁快速正反演方法及系统，来解决以上问题。
发明内容
本发明目的在于提供G-S变换的接地长导线源瞬变电磁快速正反演方法及系统，用于解决上述现有技术中存在的技术问题，利用G-S变换和OMP并行计算极大的加快了正反演技术过程，为快速获得现场接地长导线源瞬变电磁勘探数据反演结果提供了技术支撑。
为实现上述目的，本发明的技术方案是：
G-S变换的接地长导线源瞬变电磁快速正反演方法，包括以下步骤：
S1：通过汉克尔积分变换求解垂直磁场分量
Figure DEST_PATH_IMAGE001
含有贝塞尔函数的内层积分；
S2：对垂直磁场分量
Figure 508766DEST_PATH_IMAGE001
的外层采用高斯数值积分得到频率域水平地层的磁场响应；
S3：通过包含OMP并行计算的G-S算法对S2中得到的频率域水平地层的磁场响应进行频时变化，最终得到时间域的正演响应，在接地长导线源瞬变电磁正演中的频时变换部分通过G-S算法在满足精度的前提下对正演效率进行大幅度提升；
S4：在反演部分通过OMP并行计算，在Jacobi矩阵的计算中派生出若干个线程执行并行任务，最终回归迭代方程，加快反演速度。
进一步的，步骤S1具体如下：
对接地长导线源瞬变电磁法中垂直分量
Figure 997516DEST_PATH_IMAGE001
进行正演求解，首先需要通过汉克尔积分变换求解频率域垂直磁场分量
Figure 531266DEST_PATH_IMAGE001
含有贝塞尔函数的内层积分；
垂直磁场分量
Figure 789072DEST_PATH_IMAGE001
的表达式如下：
Figure 874709DEST_PATH_IMAGE002
上式中
Figure DEST_PATH_IMAGE003
为垂直磁场分量的频率域磁场响应；
Figure 319596DEST_PATH_IMAGE004
为采样角频率；
Figure DEST_PATH_IMAGE005
为电流强度；设定导线源的中心位于坐标原点，沿
Figure 860299DEST_PATH_IMAGE006
轴向两侧延伸至
Figure DEST_PATH_IMAGE007
和
Figure 159562DEST_PATH_IMAGE008
；
Figure DEST_PATH_IMAGE009
为长导线上任意一点坐标；
Figure 697991DEST_PATH_IMAGE010
为偏移距离；
Figure DEST_PATH_IMAGE011
为测点到偶极源中心点位置的距离；
Figure 958071DEST_PATH_IMAGE012
为
Figure DEST_PATH_IMAGE013
模式下的反射系数；
Figure 427099DEST_PATH_IMAGE014
为空气介质波数；
Figure DEST_PATH_IMAGE015
为线圈接收高度；
Figure 456234DEST_PATH_IMAGE016
为一阶Bessel函数；
Figure DEST_PATH_IMAGE017
为积分变量；
第一类一阶Bessel函数的汉克尔积分形式为：
Figure 431144DEST_PATH_IMAGE018
（4）
线性数字滤波公式为：
Figure DEST_PATH_IMAGE019
（5）
Figure 771995DEST_PATH_IMAGE020
（6）
上式中
Figure DEST_PATH_IMAGE021
为滤波系数点数。
进一步的，步骤S2具体如下：
在对频率域垂直磁场分量
Figure 857763DEST_PATH_IMAGE022
的内层进行积分变换求解后，对其外层通过高斯数值积分实现偶极场沿线源积分，可得到频率域水平地层的磁场响应；
高斯数值积分的形式如下：
Figure DEST_PATH_IMAGE023
（7）
将外层积分代入式（7）可得：
Figure 475826DEST_PATH_IMAGE024
（8）
上式中
Figure DEST_PATH_IMAGE025
；
Figure 74166DEST_PATH_IMAGE026
为线源长度的一半；
Figure DEST_PATH_IMAGE027
为高斯数值积节点；
Figure 715363DEST_PATH_IMAGE028
为积分系数；
Figure DEST_PATH_IMAGE029
为积分点数；采用12点高斯数值积分。
进一步的，步骤S3具体如下：
在得到频率域水平地层的磁场响应后，通过频时变换将频率域中计算得到的电磁响应转换为时间域的电磁响应；通过G-S算法进行频时变换，G-S算法是纯实数运行，需拉氏变换变量
Figure 401560DEST_PATH_IMAGE030
来替换频率域中的
Figure DEST_PATH_IMAGE031
，得到拉氏变换域中的感应电动势
Figure 326659DEST_PATH_IMAGE032
；对于给定时间
Figure DEST_PATH_IMAGE033
，感应电动势的瞬变响应
Figure 581054DEST_PATH_IMAGE034
的转换公式为：
Figure DEST_PATH_IMAGE035
（9）
Figure 37443DEST_PATH_IMAGE036
（10）
Figure DEST_PATH_IMAGE037
（11）
上式中
Figure 714281DEST_PATH_IMAGE038
为滤波系数，且
Figure DEST_PATH_IMAGE039
，
Figure 916723DEST_PATH_IMAGE040
；
Figure DEST_PATH_IMAGE041
为求和下限；
通过汉克尔积分变换、高斯数值积分、G-S频时变换算法将时间域的电磁响应推导出来，与均匀大地接地长导线源瞬变电磁线源解析公式进行对比验证；解析表达式如下：
Figure 653604DEST_PATH_IMAGE042
（12）
式（12）中，
Figure DEST_PATH_IMAGE043
为线源长度的一半；
Figure 331710DEST_PATH_IMAGE044
为自由空间磁导率常数；
Figure DEST_PATH_IMAGE045
为电导率；
Figure 766234DEST_PATH_IMAGE046
为偏移距离；
Figure DEST_PATH_IMAGE047
为测点到偶极源中心点位置的距离；
Figure 400346DEST_PATH_IMAGE048
为误差函数。
进一步的，步骤S4具体如下：
反演部分采用L1范数正则化快速反演，同时通过OMP并行计算，加快反演速度；根据正则化思想，L1范数正则化反演目标函数为：
Figure DEST_PATH_IMAGE049
（13）
式（13）中
Figure 855598DEST_PATH_IMAGE050
为数据拟合项，
Figure DEST_PATH_IMAGE051
为模型约束项，
Figure 958684DEST_PATH_IMAGE052
为正则化因子，
Figure DEST_PATH_IMAGE053
为模型向量；
数据拟合项：
Figure 915007DEST_PATH_IMAGE054
（14）
式（14）中
Figure DEST_PATH_IMAGE055
为数据加权矩阵，令其为主对角单位矩阵，
Figure 747834DEST_PATH_IMAGE056
为实测数据向量，
Figure DEST_PATH_IMAGE057
为模型正演响应函数，上标
Figure 311670DEST_PATH_IMAGE058
表示转置；
模型约束项：
Figure DEST_PATH_IMAGE059
（15）
式（15）中
Figure 89002DEST_PATH_IMAGE060
为先验模型，
Figure DEST_PATH_IMAGE061
为模型约束矩阵，
Figure 662066DEST_PATH_IMAGE062
等于
Figure DEST_PATH_IMAGE063
,下标
Figure 818241DEST_PATH_IMAGE064
表示为L1范数；
Figure DEST_PATH_IMAGE065
为最小模型约束矩阵：
Figure 271088DEST_PATH_IMAGE066
（16）
此时L1范数正则化反演目标函数可表示为：
Figure DEST_PATH_IMAGE067
（17）
由于L1范数取绝对值存在不可导的情况，取
Figure 83186DEST_PATH_IMAGE068
，将式（17）改写为：
Figure DEST_PATH_IMAGE069
（18）
将式（18）中等号右边第二项改写为：
Figure 522258DEST_PATH_IMAGE070
（19）
式（19）中
Figure DEST_PATH_IMAGE071
为对角加权矩阵，主对角元素为：
Figure 997188DEST_PATH_IMAGE072
（20）
将式（18）中模型正演响应函数
Figure DEST_PATH_IMAGE073
进行泰勒展开，取一阶线性项：
Figure 371668DEST_PATH_IMAGE074
（21）
式（21）中
Figure DEST_PATH_IMAGE075
，J为雅克比矩阵，k为迭代次数，采用中心差分方式求取：
Figure 733380DEST_PATH_IMAGE076
（22）
将式（18）关于
Figure DEST_PATH_IMAGE077
求导可得反演方程：
Figure 897514DEST_PATH_IMAGE078
（23）
式（23）可写为：
Figure DEST_PATH_IMAGE079
（24）
其中：
Figure 903647DEST_PATH_IMAGE080
（25）
Figure DEST_PATH_IMAGE081
（26）
Figure 495034DEST_PATH_IMAGE082
（27）
求解式（24）可得第k次迭代的模型修正量
Figure DEST_PATH_IMAGE083
；则第k+1次迭代计算模型为：
Figure 609621DEST_PATH_IMAGE084
（28）
正则化因子
Figure DEST_PATH_IMAGE085
取值采用以下方式：
Figure 282DEST_PATH_IMAGE086
（29）
式（29）中，
Figure DEST_PATH_IMAGE087
为第k-1次迭代的数据拟合项，
Figure 438085DEST_PATH_IMAGE088
为第k-1次迭代的模型拟合项；
设定反演结束条件为：达到最大迭代次数；相邻两次拟合差小于给定拟合误差；
Figure DEST_PATH_IMAGE089
小于给定精度
Figure 13423DEST_PATH_IMAGE090
；
拟合差计算公式为：
Figure DEST_PATH_IMAGE091
（30）。
G-S变换的接地长导线源瞬变电磁快速正反演系统，采用如上述的G-S变换的接地长导线源瞬变电磁快速正反演方法进行接地长导线源瞬变电磁快速正反演。
与现有技术相比，本发明所具有的有益效果为：
本方案其中一个有益效果在于，本发明针对地面接地长导线源瞬变电磁法，利用G-S时频变换和OMP并行运算方法相结合的方法，在正演运算过程中根据模型复杂程度将结算效率提高3-14倍，在反演运算过程中计算效率提高2倍，极大的提高了接地长导线源瞬变电磁勘探正反演计算的效率，为快速获得现场接地长导线源瞬变电磁勘探数据反演结果提供了技术支撑。。
附图说明
图1为本申请实施例的正演数值解与解析解曲线（1400m偏移距）示意图。
图2为本申请实施例的正演数值解与解析解曲线（3000m偏移距）示意图。
图3为本申请实施例的反演流程示意图。
图4为本申请实施例的OMP并行示意图。
图5为本申请实施例的多层模型反演拟合结果示意图。
具体实施方式
下面结合本发明的附图1-附图5，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
实施例：
基于G-S变换的接地长导线源瞬变电磁快速正反演方法及系统，主要包括：基于G-S变换的接地长导线源瞬变电磁快速正反演方法，以及配套的软件子系统。
基于G-S变换的接地长导线源瞬变电磁快速正反演方法流程如下：
（1）通过汉克尔积分变换求解垂直磁场分量
Figure 287410DEST_PATH_IMAGE092
含有贝塞尔函数的内层积分；
（2）对垂直磁场分量
Figure 75237DEST_PATH_IMAGE092
的外层采用高斯数值积分得到频率域水平地层的磁场响应；
（3）通过包含OMP并行计算的G-S算法对（2）中得到的频率域水平地层的磁场响应进行频时变化，最终得到时间域的正演响应，在接地长导线源瞬变电磁正演中的频时变换部分通过G-S算法在满足精度的前提下对正演效率进行了大幅度的提升；
（4）在反演部分通过单纯的OMP并行计算，在Jacobi矩阵的计算中派生出若干个线程执行并行任务，最终回归迭代方程，加快反演速度。
配套上述方法的建立基于G-S变换的接地长导线源瞬变电磁快速正反演解释软件模块包括：系统功能模块和底层支撑模块，其中，系统功能模块包括数据文件管理模块、预处理模块、快速正演成库模块、带并行反演模块和成图模块，底层支撑模块包括数据文件IO模块、嵌入式数据库模块、通用数学库模块；底层支撑模块，用于提供通用的功能函数给系统功能模块。
其中，
（1）对接地长导线源瞬变电磁法中垂直分量
Figure 836389DEST_PATH_IMAGE092
进行正演求解，首先需要通过汉克尔积分变换求解频率域垂直磁场分量
Figure 848207DEST_PATH_IMAGE092
含有贝塞尔函数的内层积分。
垂直磁场分量
Figure 875069DEST_PATH_IMAGE092
的表达式如下：
Figure DEST_PATH_IMAGE093
上式中
Figure 935429DEST_PATH_IMAGE094
为垂直磁场分量的频率域磁场响应(A/m)；
Figure DEST_PATH_IMAGE095
为采样角频率(rad/s)；
Figure 19928DEST_PATH_IMAGE096
为电流强度(A)；设定导线源的中心位于坐标原点，沿
Figure DEST_PATH_IMAGE097
轴向两侧延伸至L和-L；
Figure 874752DEST_PATH_IMAGE098
为长导线上任意一点坐标；
Figure DEST_PATH_IMAGE099
为偏移距离；
Figure 185647DEST_PATH_IMAGE100
为测点到偶极源中心点位置的距离；
Figure DEST_PATH_IMAGE101
为
Figure 236649DEST_PATH_IMAGE102
模式下的反射系数；
Figure DEST_PATH_IMAGE103
为空气介质波数；
Figure 988704DEST_PATH_IMAGE104
为线圈接收高度(m)；
Figure DEST_PATH_IMAGE105
为一阶Bessel函数；
Figure 342325DEST_PATH_IMAGE106
为积分变量。
第一类一阶Bessel函数的汉克尔积分形式为：
Figure DEST_PATH_IMAGE107
（4）
线性数字滤波公式为：
Figure 61888DEST_PATH_IMAGE108
（5）
Figure DEST_PATH_IMAGE109
（6）
上式中
Figure 401734DEST_PATH_IMAGE110
为滤波系数点数（可分为140点滤波系数和47滤波系数）；a=-3.0507817595；s=0.110599010095。
（2）在对频率域垂直磁场分量
Figure DEST_PATH_IMAGE111
的内层进行积分变换求解后，对其外层通过高斯数值积分实现偶极场沿线源积分，可得到频率域水平地层的磁场响应。
高斯数值积分的形式如下：
Figure 257563DEST_PATH_IMAGE112
（7）
将外层积分代入式（7）可得：
Figure DEST_PATH_IMAGE113
（8）
上式中
Figure 782086DEST_PATH_IMAGE114
；
Figure DEST_PATH_IMAGE115
为线源长度的一半；
Figure 739677DEST_PATH_IMAGE116
为高斯数值积节点；
Figure DEST_PATH_IMAGE117
为积分系数；
Figure 866902DEST_PATH_IMAGE118
为积分点数；采用12点高斯数值积分。12点高斯系数积分表如下：
Figure DEST_PATH_IMAGE119
（3）
在得到频率域水平地层的磁场响应后，通过频时变换将频率域中计算得到的电磁响应转换为时间域的电磁响应；本发明通过G-S算法进行频时变换，G-S算法是纯实数运行，只需极少的拉氏变换变量
Figure 921446DEST_PATH_IMAGE120
来替换频率域中的
Figure DEST_PATH_IMAGE121
，得到拉氏变换域中的感应电动势
Figure 554552DEST_PATH_IMAGE122
。对于给定时间
Figure DEST_PATH_IMAGE123
，感应电动势的瞬变响应
Figure 920812DEST_PATH_IMAGE124
的转换公式为：
Figure DEST_PATH_IMAGE125
（9）
Figure 664777DEST_PATH_IMAGE126
（10）
Figure DEST_PATH_IMAGE127
（11）
上式中
Figure 839406DEST_PATH_IMAGE128
为滤波系数
Figure DEST_PATH_IMAGE129
，选取的滤波系数
Figure 564786DEST_PATH_IMAGE130
，则G-S变换系数
Figure DEST_PATH_IMAGE131
；
Figure 496970DEST_PATH_IMAGE132
为求和下限，
Figure DEST_PATH_IMAGE133
。
通过汉克尔积分变换、高斯数值积分、G-S频时变换算法将时间域的电磁响应推导出来，与均匀大地接地长导线源瞬变电磁线源解析公式（Nabighian，1992）进行对比验证。解析表达式如下：
Figure 106942DEST_PATH_IMAGE134
（12）
式（12）中，
Figure DEST_PATH_IMAGE135
为线源长度的一半；
Figure 995133DEST_PATH_IMAGE136
为自由空间磁导率常数，其值为
Figure DEST_PATH_IMAGE137
（H/m）；
Figure 704463DEST_PATH_IMAGE138
为电导率；
Figure DEST_PATH_IMAGE139
为偏移距离；
Figure 186260DEST_PATH_IMAGE100
为测点到偶极源中心点位置的距离；
Figure 990137DEST_PATH_IMAGE140
为误差函数。
设置线源长度1.0km，发射电流1A，坐标点为（0,1400）；（0,3000），取均匀半空间为
Figure DEST_PATH_IMAGE141
，时间范围在
Figure 608200DEST_PATH_IMAGE142
等对数间隔取32个时间点，对比验证结果如图1和图2所示。
从图1和图2中可见，在同样参数情况下数值解与解析解的拟合程度好；且在
Figure DEST_PATH_IMAGE143
内相对误差在
Figure 754010DEST_PATH_IMAGE144
以内；在
Figure DEST_PATH_IMAGE145
内相对误差在
Figure 847737DEST_PATH_IMAGE146
以内；验证了G-S算法的准确性。
将G-S变换算法与传统的余弦变化算法进行单点正演时间对比，结果如表一所示：
表1 G-S算法与余弦算法单点正演时间对比
Figure DEST_PATH_IMAGE147
可以看出，随着正演模型的复杂化在满足精度要求的同时G-S算法的计算时间短，而较快的正演速度可以快速的对工区目标进行模拟，方便确定工区的测线布置，飞行高度等参数。另外，较快的正演速度会作用于反演速度，使得反演的效率提高，从而提高工作的效率，选择用G-S变换算法进行频时转换，使得正演的速度能够满足实际资料处理需要。
（4）反演部分采用L1范数正则化快速反演，同时通过OMP并行计算，加快反演速度。根据正则化思想，L1范数正则化反演目标函数为：
Figure 533933DEST_PATH_IMAGE148
（13）
式（13）中
Figure DEST_PATH_IMAGE149
为数据拟合项，
Figure 209765DEST_PATH_IMAGE150
为模型约束项，
Figure DEST_PATH_IMAGE151
为正则化因子，
Figure 713428DEST_PATH_IMAGE152
为模型向量。
数据拟合项：
Figure DEST_PATH_IMAGE153
（14）
式（14）中
Figure 904238DEST_PATH_IMAGE154
为数据加权矩阵，令其为主对角单位矩阵，
Figure DEST_PATH_IMAGE155
为实测数据向量，
Figure 331808DEST_PATH_IMAGE156
为模型正演响应函数，上标
Figure DEST_PATH_IMAGE157
表示转置。
模型约束项：
Figure 60816DEST_PATH_IMAGE158
（15）
式（15）中
Figure DEST_PATH_IMAGE159
为先验模型，
Figure 876325DEST_PATH_IMAGE160
为模型约束矩阵，
Figure DEST_PATH_IMAGE161
等于
Figure 492114DEST_PATH_IMAGE162
，下标
Figure DEST_PATH_IMAGE163
表示为L1范数。
Figure 644747DEST_PATH_IMAGE164
为最小模型约束矩阵:
Figure DEST_PATH_IMAGE165
（16）
此时L1范数正则化反演目标函数可表示为：
Figure 357488DEST_PATH_IMAGE166
（17）
由于L1范数取绝对值存在不可导的情况，可取
Figure DEST_PATH_IMAGE167
，将式（17）改写为：
Figure 16003DEST_PATH_IMAGE168
（18）
将式（18）中等号右边第二项改写为：
Figure DEST_PATH_IMAGE169
（19）
式（19）中
Figure 306038DEST_PATH_IMAGE170
为对角加权矩阵，主对角元素为：
Figure DEST_PATH_IMAGE171
（20）
将式（18）中模型正演响应函数
Figure 403307DEST_PATH_IMAGE172
进行泰勒展开，取一阶线性项：
Figure DEST_PATH_IMAGE173
（21）
式（21）中
Figure 642659DEST_PATH_IMAGE174
，
Figure DEST_PATH_IMAGE175
为雅克比矩阵，
Figure 659025DEST_PATH_IMAGE176
为迭代次数，采用中心差分方式求取：
Figure DEST_PATH_IMAGE177
（22）
将式（18）关于
Figure 249407DEST_PATH_IMAGE178
求导可得反演方程：
Figure DEST_PATH_IMAGE179
（23）
式（23）可写为：
Figure 884787DEST_PATH_IMAGE180
（24）
其中：
Figure DEST_PATH_IMAGE181
（25）
Figure 431175DEST_PATH_IMAGE182
（26）
Figure DEST_PATH_IMAGE183
（27）
求解式（24）可得第
Figure 431492DEST_PATH_IMAGE184
次迭代的模型修正量
Figure DEST_PATH_IMAGE185
。则第
Figure 305907DEST_PATH_IMAGE186
次迭代计算模型为：
Figure DEST_PATH_IMAGE187
（28）
正则化因子
Figure 931930DEST_PATH_IMAGE188
取值采用以下方式：
Figure DEST_PATH_IMAGE189
（29）
式（29）中，
Figure 83556DEST_PATH_IMAGE190
为第
Figure DEST_PATH_IMAGE191
次迭代的数据拟合项，
Figure 582671DEST_PATH_IMAGE192
为第
Figure DEST_PATH_IMAGE193
次迭代的模型拟合项。
设定反演结束条件为:①达到最大迭代次数;②相邻两次拟合差小于给定拟合误差
Figure 131333DEST_PATH_IMAGE194
;③
Figure DEST_PATH_IMAGE195
小于给定精度
Figure 718303DEST_PATH_IMAGE196
。
拟合差计算公式为:
Figure DEST_PATH_IMAGE197
（30）
反演流程图如图3所示。
OpenMP是一种共享内存系统的多处理器多线程并行语言，采用 fork-join(分叉-合并)并行执行模式。主线程作串行运算，当遇到并行模块时，调用其他从线程构成线程组，同时访问共享内存区域，执行命令，执行完毕跳出并行区域，继续执行串行命令。OpenMP并行策略具有效率高、执行快的特点，适合单机操作。接地长导线源瞬变电磁反演耗时主要在于雅克比矩阵求解时需要多次调用正演，采用二维数组存储雅可比矩阵，形成二重循环，主要运算时间也集中在内循环(每一行)参数的正演计算，参数越多耗时越长。OpenMP可以对嵌套循环体内多个循环进行并行运算，采用如图4所示方案，只对内循环进行并行运算策略，对循环内变量作无关处理，各线程对数据、函数的调用相对独立。
以H、K型地电模型为例，对使用OMP并行方法的反演程序与未使用OMP并行的反演程序进行对比，设定H型模型参数：
Figure 301600DEST_PATH_IMAGE198
，
Figure DEST_PATH_IMAGE199
，
Figure 971616DEST_PATH_IMAGE200
，
Figure DEST_PATH_IMAGE201
,
Figure DEST_PATH_IMAGE203
，
Figure 227148DEST_PATH_IMAGE204
，
Figure DEST_PATH_IMAGE205
型模型的地电参数：
Figure 929393DEST_PATH_IMAGE206
，
Figure DEST_PATH_IMAGE207
，
Figure 321191DEST_PATH_IMAGE208
，
Figure DEST_PATH_IMAGE209
,
Figure 896529DEST_PATH_IMAGE210
，
Figure DEST_PATH_IMAGE211
，初始模型均为
Figure 419783DEST_PATH_IMAGE212
的均匀半空间模型，设定模型层数为30层，每层层厚为6m，迭代次数均为40次。得到结果如下表所示：
表2 反演串并行时间对比
Figure DEST_PATH_IMAGE213
可以看出，并行时间相对串行时间效率提高了约110%，大大较少了反演迭代所需时间，实现了半航空瞬变电磁的快速反演，为半航空资料的处理节约了时间成本。
以上是本发明的较佳实施例，凡依本发明技术方案所作的改变，所产生的功能作用未超出本发明技术方案的范围时，均属于本发明的保护范围。

Claims (6)

1.G-S变换的接地长导线源瞬变电磁快速正反演方法，其特征在于，包括以下步骤：
S1：通过汉克尔积分变换求解垂直磁场分量
Figure DEST_PATH_IMAGE002
含有贝塞尔函数的内层积分；
S2：对垂直磁场分量
Figure 443569DEST_PATH_IMAGE002
的外层采用高斯数值积分得到频率域水平地层的磁场响应；
S3：通过包含OMP并行计算的G-S算法对S2中得到的频率域水平地层的磁场响应进行频时变化，最终得到时间域的正演响应，在接地长导线源瞬变电磁正演中的频时变换部分通过G-S算法在满足精度的前提下对正演效率进行大幅度提升；
S4：在反演部分通过OMP并行计算，在Jacobi矩阵的计算中派生出若干个线程执行并行任务，最终回归迭代方程，加快反演速度。
2.如权利要求1所述的G-S变换的接地长导线源瞬变电磁快速正反演方法，其特征在于，步骤S1具体如下：
对接地长导线源瞬变电磁法中垂直分量
Figure 95131DEST_PATH_IMAGE002
进行正演求解，首先需要通过汉克尔积分变换求解频率域垂直磁场分量
Figure 530660DEST_PATH_IMAGE002
含有贝塞尔函数的内层积分；
垂直磁场分量
Figure 722607DEST_PATH_IMAGE002
的表达式如下：
Figure DEST_PATH_IMAGE004
上式中
Figure DEST_PATH_IMAGE006
为垂直磁场分量的频率域磁场响应；
Figure DEST_PATH_IMAGE008
为采样角频率；
Figure DEST_PATH_IMAGE010
为电流强度；设定导线源的中心位于坐标原点，沿
Figure DEST_PATH_IMAGE012
轴向两侧延伸至
Figure DEST_PATH_IMAGE014
和
Figure DEST_PATH_IMAGE016
；
Figure DEST_PATH_IMAGE018
为长导线上任意一点坐标；
Figure DEST_PATH_IMAGE020
为偏移距离；
Figure DEST_PATH_IMAGE022
为测点到偶极源中心点位置的距离；
Figure DEST_PATH_IMAGE024
为
Figure DEST_PATH_IMAGE026
模式下的反射系数；
Figure DEST_PATH_IMAGE028
为空气介质波数；
Figure DEST_PATH_IMAGE030
为线圈接收高度；
Figure DEST_PATH_IMAGE032
为一阶Bessel函数；
Figure DEST_PATH_IMAGE034
为积分变量；
第一类一阶Bessel函数的汉克尔积分形式为：
Figure DEST_PATH_IMAGE036
（4）
线性数字滤波公式为：
Figure DEST_PATH_IMAGE038
（5）
Figure DEST_PATH_IMAGE040
（6）
上式中
Figure DEST_PATH_IMAGE042
为滤波系数点数。
3.如权利要求2所述的G-S变换的接地长导线源瞬变电磁快速正反演方法，其特征在于，步骤S2具体如下：
在对频率域垂直磁场分量
Figure DEST_PATH_IMAGE044
的内层进行积分变换求解后，对其外层通过高斯数值积分实现偶极场沿线源积分，可得到频率域水平地层的磁场响应；
高斯数值积分的形式如下：
Figure DEST_PATH_IMAGE046
（7）
将外层积分代入式（7）可得：
Figure DEST_PATH_IMAGE048
（8）
上式中
Figure DEST_PATH_IMAGE050
；
Figure DEST_PATH_IMAGE052
为线源长度的一半；
Figure DEST_PATH_IMAGE054
为高斯数值积节点；
Figure DEST_PATH_IMAGE056
为积分系数；
Figure DEST_PATH_IMAGE058
为积分点数；采用12点高斯数值积分。
4.如权利要求3所述的G-S变换的接地长导线源瞬变电磁快速正反演方法，其特征在于，步骤S3具体如下：
在得到频率域水平地层的磁场响应后，通过频时变换将频率域中计算得到的电磁响应转换为时间域的电磁响应；通过G-S算法进行频时变换，G-S算法是纯实数运行，需拉氏变换变量
Figure DEST_PATH_IMAGE060
来替换频率域中的
Figure DEST_PATH_IMAGE062
，得到拉氏变换域中的感应电动势
Figure DEST_PATH_IMAGE064
；对于给定时间
Figure DEST_PATH_IMAGE066
，感应电动势的瞬变响应
Figure DEST_PATH_IMAGE068
的转换公式为：
Figure DEST_PATH_IMAGE070
（9）
Figure DEST_PATH_IMAGE072
（10）
Figure DEST_PATH_IMAGE074
（11）
上式中
Figure DEST_PATH_IMAGE076
为滤波系数，且
Figure DEST_PATH_IMAGE078
，
Figure DEST_PATH_IMAGE080
；
Figure DEST_PATH_IMAGE082
为求和下限；
通过汉克尔积分变换、高斯数值积分、G-S频时变换算法将时间域的电磁响应推导出来，与均匀大地接地长导线源瞬变电磁线源解析公式进行对比验证；解析表达式如下：
Figure DEST_PATH_IMAGE084
（12）
式（12）中，
Figure DEST_PATH_IMAGE086
为线源长度的一半；
Figure DEST_PATH_IMAGE088
为自由空间磁导率常数；
Figure DEST_PATH_IMAGE090
为电导率；
Figure DEST_PATH_IMAGE092
为偏移距离；
Figure DEST_PATH_IMAGE094
为测点到偶极源中心点位置的距离；
Figure DEST_PATH_IMAGE096
为误差函数。
5.如权利要求4所述的G-S变换的接地长导线源瞬变电磁快速正反演方法，其特征在于，步骤S4具体如下：
反演部分采用L1范数正则化快速反演，同时通过OMP并行计算，加快反演速度；根据正则化思想，L1范数正则化反演目标函数为：
Figure DEST_PATH_IMAGE098
（13）
式（13）中
Figure DEST_PATH_IMAGE100
为数据拟合项，
Figure DEST_PATH_IMAGE102
为模型约束项，
Figure DEST_PATH_IMAGE104
为正则化因子，
Figure DEST_PATH_IMAGE106
为模型向量；
数据拟合项：
Figure DEST_PATH_IMAGE108
（14）
式（14）中
Figure DEST_PATH_IMAGE110
为数据加权矩阵，令其为主对角单位矩阵，
Figure DEST_PATH_IMAGE112
为实测数据向量，
Figure DEST_PATH_IMAGE114
为模型正演响应函数，上标
Figure DEST_PATH_IMAGE116
表示转置；
模型约束项：
Figure DEST_PATH_IMAGE118
（15）
式（15）中
Figure DEST_PATH_IMAGE120
为先验模型，
Figure DEST_PATH_IMAGE122
为模型约束矩阵，
Figure DEST_PATH_IMAGE124
等于
Figure DEST_PATH_IMAGE126
,下标
Figure DEST_PATH_IMAGE128
表示为L1范数；
Figure DEST_PATH_IMAGE130
为最小模型约束矩阵：
Figure DEST_PATH_IMAGE132
（16）
此时L1范数正则化反演目标函数可表示为：
Figure DEST_PATH_IMAGE134
（17）
由于L1范数取绝对值存在不可导的情况，取
Figure DEST_PATH_IMAGE136
，将式（17）改写为：
Figure DEST_PATH_IMAGE138
（18）
将式（18）中等号右边第二项改写为：
Figure DEST_PATH_IMAGE140
（19）
式（19）中
Figure DEST_PATH_IMAGE142
为对角加权矩阵，主对角元素为：
Figure DEST_PATH_IMAGE144
（20）
将式（18）中模型正演响应函数
Figure DEST_PATH_IMAGE146
进行泰勒展开，取一阶线性项：
Figure DEST_PATH_IMAGE148
（21）
式（21）中
Figure DEST_PATH_IMAGE150
，J为雅克比矩阵，k为迭代次数，采用中心差分方式求取：
Figure DEST_PATH_IMAGE152
（22）
将式（18）关于
Figure DEST_PATH_IMAGE154
求导可得反演方程：
Figure DEST_PATH_IMAGE156
（23）
式（23）可写为：
Figure DEST_PATH_IMAGE158
（24）
其中：
Figure DEST_PATH_IMAGE160
（25）
Figure DEST_PATH_IMAGE162
（26）
Figure DEST_PATH_IMAGE164
（27）
求解式（24）可得第k次迭代的模型修正量
Figure DEST_PATH_IMAGE166
；则第k+1次迭代计算模型为：
Figure DEST_PATH_IMAGE168
（28）
正则化因子
Figure DEST_PATH_IMAGE170
取值采用以下方式：
Figure DEST_PATH_IMAGE172
（29）
式（29）中，
Figure DEST_PATH_IMAGE174
为第k-1次迭代的数据拟合项，
Figure DEST_PATH_IMAGE176
为第k-1次迭代的模型拟合项；
设定反演结束条件为：达到最大迭代次数；相邻两次拟合差小于给定拟合误差；
Figure DEST_PATH_IMAGE178
小于给定精度
Figure DEST_PATH_IMAGE180
；
拟合差计算公式为：
Figure DEST_PATH_IMAGE182
（30）。
6.G-S变换的接地长导线源瞬变电磁快速正反演系统，其特征在于，采用如权利要求1-5任一项所述的G-S变换的接地长导线源瞬变电磁快速正反演方法进行接地长导线源瞬变电磁快速正反演。
CN202211243244.3A 2022-10-11 2022-10-11 G-s变换的接地长导线源瞬变电磁快速正反演方法及系统 Active CN115508898B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202211243244.3A CN115508898B (zh) 2022-10-11 2022-10-11 G-s变换的接地长导线源瞬变电磁快速正反演方法及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202211243244.3A CN115508898B (zh) 2022-10-11 2022-10-11 G-s变换的接地长导线源瞬变电磁快速正反演方法及系统

Publications (2)

Publication Number Publication Date
CN115508898A true CN115508898A (zh) 2022-12-23
CN115508898B CN115508898B (zh) 2023-11-10

Family

ID=84510845

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202211243244.3A Active CN115508898B (zh) 2022-10-11 2022-10-11 G-s变换的接地长导线源瞬变电磁快速正反演方法及系统

Country Status (1)

Country Link
CN (1) CN115508898B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116879964A (zh) * 2023-08-14 2023-10-13 成都理工大学 一种时频电磁频率域数据自约束稳健电阻率反演方法

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102419453A (zh) * 2011-07-15 2012-04-18 中国科学院地质与地球物理研究所 长导线源瞬变电磁地空探测方法
CN104237956A (zh) * 2014-03-06 2014-12-24 长安大学 电性源瞬变电磁地空探测方法
US20200003928A1 (en) * 2018-07-02 2020-01-02 Institute Of Geology And Geophysics, Chinese Academy Of Sciences Method and device for surface-borehole transient electromagnetic detection
US20200025958A1 (en) * 2017-03-29 2020-01-23 Westerngeco Llc Compressive sensing imaging
CN111965714A (zh) * 2020-07-15 2020-11-20 中国地质大学（武汉） 一种基于暂态过程的电磁探测方法、设备及存储设备
CN113933905A (zh) * 2021-09-30 2022-01-14 中国矿业大学 一种圆锥型场源瞬变电磁反演方法
CN215894978U (zh) * 2021-08-05 2022-02-22 成都理工大学 基于g-s变换的半航空瞬变电磁成像系统
CN115047530A (zh) * 2022-06-15 2022-09-13 吉林建筑大学 基于一维反演的地空瞬变电磁数据三维频率域解释方法

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102419453A (zh) * 2011-07-15 2012-04-18 中国科学院地质与地球物理研究所 长导线源瞬变电磁地空探测方法
CN104237956A (zh) * 2014-03-06 2014-12-24 长安大学 电性源瞬变电磁地空探测方法
US20200025958A1 (en) * 2017-03-29 2020-01-23 Westerngeco Llc Compressive sensing imaging
US20200003928A1 (en) * 2018-07-02 2020-01-02 Institute Of Geology And Geophysics, Chinese Academy Of Sciences Method and device for surface-borehole transient electromagnetic detection
CN111965714A (zh) * 2020-07-15 2020-11-20 中国地质大学（武汉） 一种基于暂态过程的电磁探测方法、设备及存储设备
CN215894978U (zh) * 2021-08-05 2022-02-22 成都理工大学 基于g-s变换的半航空瞬变电磁成像系统
CN113933905A (zh) * 2021-09-30 2022-01-14 中国矿业大学 一种圆锥型场源瞬变电磁反演方法
CN115047530A (zh) * 2022-06-15 2022-09-13 吉林建筑大学 基于一维反演的地空瞬变电磁数据三维频率域解释方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
HAIYAN YANG ET AL.: "Transient Electromagnetic Response With a Ramp Current Excitation Using Conical Source", 《IEEE》, vol. 7, pages 63829 - 63836, XP011726886, DOI: 10.1109/ACCESS.2019.2914740 *
杜兴忠 等: "航空瞬变电磁正演中频时转换方法分析及系数的选取", 《物探化探计算技术》, vol. 37, no. 6, pages 671 - 679 *
罗延钟 等: "G-S变换的快速算法", 《地球物理学报》, vol. 43, no. 5, pages 684 - 690 *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116879964A (zh) * 2023-08-14 2023-10-13 成都理工大学 一种时频电磁频率域数据自约束稳健电阻率反演方法
CN116879964B (zh) * 2023-08-14 2024-04-26 成都理工大学 一种时频电磁频率域数据自约束稳健电阻率反演方法

Also Published As

Publication number Publication date
CN115508898B (zh) 2023-11-10

Similar Documents

Publication Publication Date Title
US11460519B2 (en) Method for making a magnetic gradiometer with high detection accuracy and success rate
Sailhac et al. Identification of sources of potential fields with the continuous wavelet transform: Two‐dimensional wavelets and multipolar approximations
CN107632964B (zh) 一种平面地磁异常场向下延拓递归余弦变换法
CN111399066B (zh) 一种基于正交基函数处理标量磁异常梯度信号的方法
CN110133644B (zh) 基于插值尺度函数法的探地雷达三维正演方法
Wang et al. Seismic velocity inversion transformer
CN110687610A (zh) 一种基于重磁数据相关性分析的场源定位及属性识别方法
CN115508898A (zh) G-s变换的接地长导线源瞬变电磁快速正反演方法及系统
Zhang et al. 3-D image-domain least-squares reverse time migration with L1 norm constraint and total variation regularization
Lou et al. Seismic volumetric dip estimation via multichannel deep learning model
Wang et al. 2D joint inversion of CSAMT and magnetic data based on cross-gradient theory
Gao et al. An efficient vector elastic reverse time migration method in the hybrid time and frequency domain for anisotropic media
CN107942374A (zh) 绕射波场提取方法和装置
Tobely et al. Position detection of unexploded ordnance from airborne magnetic anomaly data using 3-D self organized feature map
CN116226627B (zh) 一种非高斯环境下洛伦茨约束角度估计方法及系统
Chung et al. An adaptive phase space method with application to reflection traveltime tomography
CN103513288B (zh) 一种二维网格数据的补偿方向滤波方法
CN109387872A (zh) 表面多次波预测方法
Zhang et al. Preconditioned transmission+ reflection joint traveltime tomography with adjoint-state method for subsurface velocity model building
Yang et al. Determination of the phase-velocity directions in anisotropic media using a direction vector
Tang et al. Subsalt velocity estimation by target-oriented wave-equation migration velocity analysis: A 3D field-data example
Singer et al. Modelling of electromagnetic fields in thin heterogeneous layers with application to field generation by volcanoes—theory and example
Ravve et al. Eigenrays in 3D heterogeneous anisotropic media: Part IV--Geometric spreading from traveltime Hessian
Bernasconi et al. 3-D traveltimes and amplitudes by gridded rays
Liu et al. Research on a matching detection method for magnetic anomaly of underwater targets

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