基于Matlab平台的地-井TEM快速成像

杨 毅， 张 杰， 霍美净， 邓晓红，王兴春， 智庆全， 武军杰

(1.中国地质科学院 地球物理地球化学勘查研究所，廊坊 065000;2．国土资源部 地球物理电磁法探测技术重点实验室，廊坊 065000;3．成都理工大学，成都 610059;4.廊坊市水利勘察规划设计院，廊坊 065000)

0 引言

地-井TEM法是瞬变电磁法的一个重要分支方法，兴起于上世纪70年代，最早由苏联、加拿大等国开展相关研究。Crone于上世纪70年代研制了第一套井中瞬变电磁测量系统，标志着这项方法技术应用于实际勘察工作。多年来经过国内外学者的努力，不论是仪器系统[1-3]还是方法技术[4-12]都得到了极大发展，使得地-井TEM法广泛应用于国内外的硫化物型铜镍矿的勘察[13-23]，并取得巨大成功。

1 理论方法

1.1 等效涡流理论

地-井TEM法一般采用地面发射，井中接收的工作方式，能以“一孔之见”对井旁和井底可能存在的异常做出准确预判，为进一步的地质钻孔设计提供重要参考。井中瞬变电磁法最早是通过研究自由空间中局部异常体(薄板体、球体等)的特征来开展的[24-25]，其场的特征研究均以感应涡流为基础，其中又以导电薄板(图1)最具代表性，导电薄板内感应涡流建立与消失的物理过程，为等效电流环反演奠定理论基础。

图1 导电薄板内涡旋电流分布示意图Fig.1 Schematic diagram of eddy current distribution in conductive thin plate

图2 圆形载流环参数模型Fig.2 Parametric model of circular current

位于一次场中的导电薄板，在发射突然关断瞬间，根据法拉第定律，为了维持关断之前一次场的作用，在导电薄板体内会产生感应涡流以维持关断之前一次场的作用，感应涡流在板内产生与薄板体形态相近的一系列电流强度和大小均不同的电流环，对应于瞬变过程早期的电流环主要集中在板体的边缘，尺寸相对大，这些电流环随着时间推移逐渐向导体中心扩散，直至以热损耗形式消亡。在这过程中导电薄板体内的任意时刻的感应电流分布可以用一个等效电流环表示，如图2所示。这样通过等效电流环即实现了对局部导电薄板的描述，位于自由空间的电流环可由七个参数来表示，即电流环中心坐标(x，y，z)、电流环半径、电流强度、电流环倾角、电流环方位角。

等效载流环包含七个参数(图2)，其在自由空间中任意点磁场的数学表达式为：

(1)

(2)

依据式 ( 1 ) 和式 ( 2 ) 即可计算自由空间中任意电流强度、大小、展布的电流环在任一点产生的场。

1.2 异常场提取

在异常场成像之初需要对实测数据进行异常提取，采用层状大地响应来近似等效背景围岩的响应，从实测数据中减去这个背景场，从而得到异常场响应。计算上采用Kerry Key(2009)[26]的方法，从电偶极子出发推导其在层状大地空间中任意点的电磁响应，再通过有限个长度的偶极子叠加获得回线源层状大地瞬变电磁响应(李建平，2005)[27]。

水平电偶极子位于层状大地时，除z方向的矢量位之外，还出现x方向的分量，其柱坐标下矢量势分量分别写为式(3)与式(4)。

(3)

(4)

其相应的磁场各分量表达式可通过式(5)～式(7)求得。

(5)

(6)

(7)

通过以上步骤求得的磁场响应为包含贝塞尔函数的复杂积分，求解上采用Hankel数字滤波计算，其一般表达式可写为式(8)。

(8)

完成频率域场的计算再采用余弦变换将场值换算到时间域，以Hz为例其数值离散近似表达式写为式(9)。

(9)

依据以上公式，使用Fortran编写层状大地正演计算代码，为人机交互功能控件提供回调函数。

2 人机交互快速成像实现

Matlab图形用户界面 ( GUI ) 是指由窗口、菜单、图标、光标、按键、对话框和文本等各种图形对象组成的用户界面。它让用户定制用户与Matlab的交互方式。用户通过鼠标或键盘选择、激活这些图形对象，使计算机产生某种动作或变化。基本图形对象分为控件对象和用户界面菜单对象，简称控件和菜单。

Matlab提供了一套可视化的创建图形窗口的工具，使用图形用户界面开发环境可方便地创建GUI应用程序，它可以根据用户设计的GUI布局，自动生成M文件的框架，用户使用这一框架编制自己的应用程序。Matlab提供了一套可视化的创建图形用户接口 ( GUI ) 的工具，层次架构如图3。

图3 GUI功能及程序运行流程简图Fig.3 GUI function and program flow diagram

2.1 人机交互界面设计和功能实现

交互界面的布局应该在保证基本功能的前提下尽量简洁美观，秉承这一理念，根据需要将快速成像程序分成三个功能模块来实现①数据载入编辑；②人机交互背景场拟合；③人机交互异常场拟合。

其功能简图及其主代码名称如下图3所示，运行主程序PureAnomalyFitting.m文件，点击菜单Data Raw Process调用inversion.m程序，进入界面，在此界面完成三分量地-井TEM数据的读入和编辑工作，完成后关闭子GUI界面将数据传递给主GUI。主GUI界面左侧分区可进行地层参数、发射框角点坐标的输入、发射电流输入、下降沿输入以及场值计算类型的选择，输入完毕后，绘制地层厚度和电阻率折线图，程序中的交互设计允许在此图上对层厚和电阻率值进行拖动，快速实现层厚和电阻率参数值修改，再进行层状大地响应正演拟合，最后保存背景场和异常场。主GUI界面右侧分区可以进行电流环个数以及七参数的快速载入，点击电流环正演，进行异常场拟合，同样在另一个图形中显示电流环空间位置展布图，在此设计了交互编辑功能，允许用户在不同视角下对电流环的参数进行快速编辑，然后根据异常特征进行电流环交互拟合，人工试错不断修正，直至达到最佳拟合，最后将异常场不同道拟合最好的电流环在三维坐标里成图即可得到异常的三维空间展布。

2.2 数据载入编辑

程序开始部分首先应该对实测数据进行载入并以剖面曲线的形式展现出来，为了美观和数据操作编辑方便，将其作为一个子GUI单独呈现，其作用为为初始成像做好基本数据处理，包含数据载入、数据编辑、曲线圆滑等三个大的功能。代码主M文件为inversion.m，对应的GUI文件名为inversion.fig。

GUI运行后界面如图4所示，设计了三个axes控件，分别用以显示地-井TEM测量的x、y、z三个分量数据，在Data Raw Process下拉菜单PLOT选项可对需要显示的道数进行选择。菜单的安排上，设置了Data Load、Data Plot、Data Smooth、Data Delete四个功能，分别用以完成数据载入、成图、数据圆滑、数据点删除等功能，其中在plot菜单下调用selectchannel子GUI用于选择需要操作的时间道。

图4 子GUI运行效果示意图Fig.4 Sketch map of running effect of sub GUI

图5 主GUI控件设计示意图Fig.5 Schematic diagram of the main GUI control

2.3 背景场拟合和异常场提取

M代码主文件名PureAnomalyFitting.m文件，对应的GUI文件文PureAnomalyFitting.fig，其布设如图5所示，按功能分区，左侧为背景场拟合部分，使用的控件有一个Button Group，内设两个互斥的Radio Button，分别代表背景场计算场值为电动势还是磁场；Layer earth H&R负责层状大地层厚和层电阻率的录入，具体方式为在界面上手动输入；TXI ( A )和RAMP(ms)分别对应发射电流和下降沿；Loop corner coordinates负责发射框角点坐标的录入，可在界面上手动输入；Axes6负责显示录入地层的厚度和电阻率(折线图)；Axes1、Axes2、Axes3分别用于拟合曲线显示，axes4用于拟合差的显示，这几个axes控件为公共控件，在异常场拟合部分也用于拟合曲线和拟合差的显示。pure anomaly按钮用于拟合曲线和拟合差曲线的绘制。

对层状大地层厚和电阻率的交互编辑功能为本段程序的重点实现部分，在初始给定层数、层厚度和电阻率之后，程序允许在Axes6界面中对顶底界面、电阻率折线图进行拖动，实现层状大地层厚和电阻率参数的快速修改。

在菜单中设置integral菜单，其功能为在完成层状大地正演后，从实测数据中减去层状大地响应，获取异常场响应，并对异常场积分，获取磁场响应，为异常场的电流环成像做好准备。

2.4 异常场拟合

在完成背景场拟合和异常场提取之后即进行此步骤，在GUI布局设计上，如图5所示，从右至左，filament choose为下拉菜单控件，用于电流环参数设置后各参数的载入；下拉菜单之下为七组两两横排的控件，用于电流环七参数的输入，左侧为edit text可编辑文本框控件，右侧为slide滑块控件，滑块和文本框输入设置为联动，即拖动滑块或输入文本框，另一项值也随着同步变化，可以快速实现参数的修改；Edit Panel Select 控件是一个互斥的radio button组，用于交互操作时，交互区域的选择，分别对应背景场拟合时地层层厚、电阻率折线图和电流环两个分区的选择和锁定。Axes5用于发射框、井、一次场以及给定电流环的空间位置三维显示；filament plot按钮控件用于初次电流环的绘制以及参数重置；refresh按钮用于交互成像电流环参数改变后对电流环的刷新显示。

对电流环的交互编辑功能为本段程序的重点实现部分，在初始给定电流环个数和基本位置之后，程序允许在x-y，x-z，y-x三个视角界面下对电流环位置进行拖动，实现参数的快速修改。

2.5 GUI主程序运行流程及其功能实现

本程序的编制基于Matlab的GUI平台，根据需要设置了控件，并逐一完善控件功能。下面将简单叙述笔者编写的人机交互程序的工作流程，实现功能的控件、菜单以及回调函数。

打开Matlab程序，将路径切换到目标文件夹，双击PureAnomalyFitting.m函数，点击运行按钮主程序开始运行，弹出如下界面(图6)。

1 ) 完成数据的读入与编辑。单击Data Raw Process菜单，弹出inversion.fig界面(图7)，在此界面下载入数据，并可对数据进行简单编辑，其调用的函数文件为loaddata1.m，主要完成实测数据文件的整理和读入；Data Raw Process点击后调用inversion.fig负责数据的初步编辑处理，其调用的函数文件为loaddata1.m，其主要功能为实测数据文件的整理和读入。

2 ) 完成发射参数和层状大地参数给定。从Button Group中选择层状大地正演计算(背景场)所需场类型，df为感应电动势，ff为感应磁场；在TXI和RAMP两个可编辑输入文本控件中分别输入发射电流和下降沿；在layer earth H&R输入文本框中输入正演电阻率模型参数；在Loop corner coordinates可编辑文本框里输入发射框角点坐标；完成之后点击Layer earth model plot菜单，即可在GUI左下角的axes6中绘制出电阻率-层厚折线图；点击Primary field plot可在axes5中绘制一次场矢量图，发射框和井的位置图，主要调用ZY.m完成一次矢量场的空间分布计算(图8)。

3 )通过人机交互方式完成背景场计算。根据实测曲线特征，已经完成了层状大地正演所需参数准备，点击Layer earth forward按钮完成层状大地背景场正演计算，此处主要调用temFwdLayerBoreholeFast.exe执行文件；在完成正演后，按下pure anomaly按钮绘制实测数据与层状大地正演拟合曲线(图9)；如果拟合不满意，则可以点GUI击右下角Edit panel select下的layer earth edit单选按钮，启动GUI左下角axes6上的电阻率-层厚折线图编辑模式，在axes6范围内左键选中线段(图9)，线段被选中后即可进行快速拖动以改变电阻率或层厚度值(水平方向线段可上下拖动，改变地层厚度，垂直方向线段可左右拖动，改变层电阻率值)，然后再点击 Layer earth forward按钮进行正演，直到拟合满意为止。

4 ) 完成异常场转换。点击integral菜单，用“差值法”完成异常场提取并将电动势积分转换为磁场响应，调用函数jifen11.m。

图6 GUI主程序运行界面Fig.6 GUI main program running interface

图7 副GUI程序运行界面Fig.7 Vice GUI program running interface

图8 主GUI程序分区参数输入示意图Fig.8 Schematic diagram of main GUI program partition parameter input

图9 主GUI程序背景拟合示及地层参数编辑意图Fig.9 Background fitting of main GUI program and editing intention of formation parameters

图10 异常场拟合与电流环编辑示意图Fig.10 Abnormal field fitting and current loop edit diagram

5 ) 完成电流环初始参数值的给定。点击filamentNumber菜单，在弹出对话框中输入参与异常场拟合的电流环个数，点击后载入到下拉菜单filament choose中，然后在其正下方表征电流环七参数的七个可编辑输入文本框中给出每个电流环的参数，注意每完成一个电流环参数的输入，需要在下拉菜单中选中对应的电流环名称已完成参数的载入，完成后点击Filament plot按钮将电流环绘制于axes5中完成显示。

6 ) 通过人机交互方式完成电流环拟合。点击GUI右下角的Filament Edit radio控件选中电流环编辑模式，则可在x-y，x-z，y-z三个平面视角内对电流环进行拖动操作，完成其空间位置坐标的迅速修改，三个平面视角的实现由三个activity 控件控制，对电流环大小、强度以及两个空间角度的修改则可以在图中点击选中(图10)相应的电流环，然后在电流环输入文本里编辑或拖动其后的滑块实现，完成后点击GUI下方的refresh按钮完成数据更新以及电流环的重绘。

完成电流环设置后，可点击pure plot按钮，在axes1、axes2、axes3、axes4中分别显示各分量异常场与电流环响应拟合情况(图10)，如不满意则可返回上一步采用人机交互方式完成电流环参数的快速时时修正。

7 ) 数据的输出与显示。在完成各道异常场曲线与对应电流环的拟合后，点击spatial按钮，将不同道拟合得到的电流环会绘制于一个图上，最后用800 dpi分辨率将图从四个视角以*.jpeg格式输出。至此，完成了整个交互成像程序的编制。

3 理论算例

使用Maxwell软件中设计的理论模型计算结果作为假设实测数据进行交互成像。

模型均匀半空间中放置一个薄板体，其中均匀半空间电阻率为300 Ω·m，中心框发射，框边长为200 m，采样道为36道，间隔采用Crone 50 ms间隔，下降沿为0.5 ms，发射电流10 A，计算单位为pT/s；导电板体大小设置为50 m×50 m，纵向电导为500 S，导电薄板中心坐标为(-50,-50，-100)，倾角为25°，旋转角为30°。

图11 异常提取与电流环成像拟合曲线Fig.11 Fitting curve of anomaly extraction and current loop imaging

采用上一节所述方法对Maxwell设计的半空间含导电薄板模型计算的结果进行交互成像，其中背景场拟合(图11(a))在早期相对差，在中晚期拟合较好，自编程序计算的半空间响应能够完全反映背景场信息，通过“差值法”从实测数据中减去自编程序模拟计算的背景场即得到异常场，并采用电流环响应与之拟合(图11(b))，对于三个分量拟合结果在早期道相对差，自编程序在给定与Maxwell背景值一致条件下计算得到的场值相对强，原因可能为未考虑背景与导电薄板的耦合所致，中晚期拟合好(图11(b)、图11(c)、图11(d))，所得电流环与薄板大小相近，形态一致(图12)。表明了程序的正确性，为下一步实测数据的交互成像奠定了基础。

4 结论

笔者实现了基于等效电流环理论的Matlab 人机交互快速成像，编写了人机交互程序，主要取得如下结论：

1)编写了基于Matlab GUI平台的等效电流环快速成像界面，利用GUI自带的动作控件函数实现了层状大地模型和电流环模型的快速修改编辑，为异常场提取和成像奠定了基础。

2)完成了理论数据的快速成像，所得电流环与理论导电薄板模型一致，验证了程序的正确性，为地-井TEM实测数据的定量解释开辟了新思路。

图12 等效电流环空间分布图Fig.12 Spatial distribution diagram of equivalent current