EH4数据处理在SCS2D二维反演平台中的实现

李美英+屈栓柱

摘 要 Matlab高级语言开发平台的函数资源丰富，对矩阵的计算简便快捷，文件输出方便。以这些特点为契机，在Matlab平台上，转换EH4仪器采集的原始Z文件为SCS2D二维反演所需的MTD文件，可大大提高EH4数据的处理能力。

关键词 音频大地电磁法；Matlab平台；SCS2D反演平台；代码

中图分类号：TP3 文献标识码：A 文章编号：1671-7597（2014）05-0041-03

可控源音频大地电磁法（CSAMT）是在大地电磁法（MT）和音频大地电磁法（AMT）的基础上发展起来的一种人工源频率测深方法。大地电磁法（MT）的观测频率很低（n×10-5～n×104 Hz），所以勘探深度很大，可达100千米以上。但由于频率偏低，所以对浅层的分辨率较差，而且生产效率也比较低。为了更好地研究几十米到几千米深度范围内的地电结构，在MT方法的基础上，形成了音频大地电磁法（AMT）。音频大地电磁法的工作方法、观测参数与MT方法相同，不过该方法主要观测由于雷电作用产生的音频（n×100～n×104 Hz）大地电磁场，由于音频大地电磁法的观测频段内，天然大地电磁场的强度较弱，人文干扰大，所以信噪比很低，野外数据采集十分困难，需要长时间采集并采用多次叠加技术。为了克服上述困难，20世纪70年代初，加拿大多伦多大学的D.W.Strangway教授和他的学生提出沿用AMT的观测方式，观测人工供电产生的音频电磁场。由于所观测电磁场的频率、强度和方向可以由人工控制，而其观测方式又与AMT法相同，所以叫可控源音频大地电磁法。

EH4连续电导率剖面仪是由美国EMI公司生产的电磁法仪器系统，属于混合源方法的仪器。它利用的场源可以是天然场，也可以是人工场。数据采集方式和AMT法相同；在频率（10～1000）Hz范围内，采集天然电磁场信号；在频率（750～92K）Hz范围内，天然场高频成分信号比较弱，使用人工源信号，这时满足波区条件的收发距较小，容易实现；并且发送信号的时间短，功率小，装备和电源轻便。本仪器配备高频磁探头，观测的电磁信号频率范围为（10～92k）Hz，可用于测量地下几米至1千米范围内介质电阻率的连续分布情况。

近十年来，国际上出现了一批代表性的电磁法数据处理解释软件，如美国Zonge公司开发的SCS2D反演软件、意大利Geosystem公司开发的一个综合性解释平台WinGlink、吉林大学开发的GeoElectro电法数据处理系统和美国某大学实验室研发的Sinv2d反演软件等，不同反演软件产生的反演模型拟合差因数据量及侧重深度大小不一，使用一种成熟的处理软件不仅便于资料对比，同时也为资料的地质解译提供了较好的物探参考。

随着电磁法的普遍应用，加拿大凤凰公司的V8仪器越来越被新疆众多地质单位青睐，配套的SCS2D反演软件也逐渐展露其头角，本文利用第四代计算机语言Matlab为平台，对EH4仪器采集的数据格式进行转换，使之能够通过SCS2D反演，大大提高了数据处理的能力。

1 数据格式

软件所涉及的文件包括EH4仪器自身采集并初步处理取得的Z文件（文件个数与点位相同）、数据处理所需的参数文件*.par、生成供SCS2D反演的*.mtd文件。

1）参数文件。参数文件（*.par）实际上是一个索引目录文件，它展示了实际点位、高程与该点位所测原始数据的对应关系。文件共五列，分别为点号、坐标x、坐标y、高程、文件号，可用文本进行编辑。

2）Z文件格式。EH4仪器首先采集的是人工和天然大地电磁场强度变化的时间序列，然后内部进行傅立叶变换得出频率域幅值响应再通过卡尼亚电阻率公式（见式（1））计算后把结果写入Z文件，其内容包括频点、ExHy的相干系数、ExHy标量视电阻率、ExHy标量相位、EyHx的相干系数、EyHx标量视电阻率、EyHx标量相位及各方向阻抗张量的实部及虚部等。

（1）

式中：f—频率（Hz）；ρ-视电阻率（Ω·m）。

3）*.MTD文件格式。*.MTD是测线所有点位参数的集成文件，用逗号分隔，遗漏数据在行中用双逗号表示。如图1，由Z文件转换而来的数据主要覆盖5～9列，1～4列根据参数文件（*.par）直接写入。1～9列表示的含义见表1。

2 Matlab平台下的数据格式转换过程

Matlab是矩阵实验室（Matrix Laboratory）的简称，是美国MathWorks公司出品的商业数学软件，用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境。

1）友好的工作平台和编程环境。Matlab由一系列工具组成。这些工具方便用户使用Matlab的函数和文件，其中许多工具采用的是图形用户界面（图2）。包括Matlab桌面和命令窗口、历史命令窗口、编辑器和调试器、路径搜索和用于用户浏览帮助、工作空间、文件浏览器等。随着Matlab的商业化以及软件本身的不断升级，Matlab的用户界面也越来越精致，更加接近Windows的标准界面，人机交互性更强，操作更简单。而且新版本的Matlab提供了完整的联机查询、帮助系统，极大的方便了用户的使用。简单的编程环境提供了比较完备的调试系统，程序不必经过编译就可以直接运行，而且能够及时地报告出现的错误及进行出错原因分析。

图2 程序运行界面

2）简单通用的程序语言。Matlab使用的第四代计算机语言包含控制语句、函数、数据结构、输入和输出语句等，具有面向对象编程的特点。用户可以在命令窗口中将输入语句与执行命令同步，也可以先编写好一个较大的复杂的M文件后再一起运行（图3为已编好的M文件运行界面）。新版本的Matlab语言是基于最为流行的C++语言基础上的，因此语法特征与C++语言极为相似，而且更加简单，更加符合科技人员对数学表达式的书写格式。使之更利于非计算机专业的科技人员使用。而且这种语言可移植性好、可拓展性极强，这也是Matlab能够深入到科学研究及工程计算各个领域的重要原因。endprint

图3 运行程序调用参数文件

3）软件流程。野外EH4仪器数据采集完毕，测量台班会记录并形成一个关于点号、坐标、高程、文件号的索引参数文件，根据该文件读取原始z文件中xy方向的相关系数、视电阻率、相位、频率等，Matlab程序会根据提供的各种数据组合成SCS2D反演软件所需的mtd格式文件，并最终输出。软件大致流程如图4。

图4 软件流程图

4）程序代码及注释。

clear all；clc；

[cfilename，cdir]=uigetfile（{'*.par'，'parameter file（*.par）'；'*.*'，'All files（*.*）'}，...

'Choose a parameter file to read'）；%选择参数文件并读取。

cfile=[cdir，cfilename]；

fid=fopen（cfile，'r'）；

data=fscanf（fid，'%f'，[5，inf]）；%参数文件分四列，为点号，坐标x，坐标y，高程，文件号。

data=data'；

nsite=size（data，1）；

j=1；

for isite=1：nsite%进入循环语句。

sname=num2str（data（isite，5））；

slength=length（sname）；

if slength==1

sname=['00'，sname]；

elseif slength==2

sname=['0'，sname]；%确定z文件名。

end

zfile=['zhs.'，sname]；

zfid=fopen（zfile，'r'）；

zdata=fscanf（zfid，'%f，[15，inf]）；%读取z文件

zdata=zdata'；

for i=1：39；

if zdata（i，3）<=0；

continue；

else

mdata（j，6）=zdata（i，3）；

mdata（j，5）=zdata（i，1）；

mdata（j，1：4）=data（isite，1：4）；

mdata（j，7）=0.5；

mdata（j，8）=zdata（i，4）*3.1415926*1000/180；

mdata（j，9）=0.5；

mdata（j，10：13）=0；

mdata（j，14）=mdata（j，6）；

mdata（j，15）=mdata（j，8）；

mdata（j，16：17）=0；%生成mtd文件矩阵

j=j+1；

end

end

fclose（zfid）；

end

fclose （fid）；

mtdfile=['line'，num2str（data（1，2）），'.mtd']；

fid=fopen（mtdfile，'w'）；

fprintf（fid，'%s＼n'，'"Stn"，"GridE"，"GridN"，"Elev"，"Freq"，"ARTMobs"，"ARTMerr"，"ZPTMobs"，"ZPTMerr"，"ARTEobs"，"ARTEerr"，"ZPTEobs"，"ZPTEerr"'）；

fprintf（fid，'%s＼n'，'"From S2D2D v2.10l Date：16/12/09 Time：13：34：14"'）；

for i=1：j-1

fprintf（fid，'%6.1f，%6.1f，%6.1f，%6.1f，%6.1f，%6.1f，%6.1f，%6.1f，%6.1f'，mdata（i，1：9））；

fprintf（fid，'%s'，'，，，，'）；

fprintf（fid，'%6.1f，%6.1f'，mdata（i，14：15））；

fprintf（fid，'%s＼n'，'，，，'）；%输出mtd文件

end

fclose（fid）；

3 实例分析

1）工区概况。库普边防派出所位于奇台县八一牧场西北侧，距奇台县城约197千米，为昌吉边防支队的下属单位，多年来一直处于缺水状态。为做好双拥工作，更为解决官兵的实际用水困难，应昌吉国土资源局邀请，物探队派遣EH4项目组，了解该区地下电阻率分布特征，为寻找储水构造并布置钻井提供依据。

该区地处东准噶尔碱性花岗岩体中，寻找裂隙水是本次工作的主攻目标。

2）工作布置。如图5，限于工区实际条件，布置测线：3线、4线、5线。取点数较多的3线为例，3线307点西北侧为学校办公楼，楼前308点西约6米处为已知水井，该水井为前些年所打，井深10米，基底为花岗岩。由于上游水量降低，该井出水量日益减少，已不能满足学校最基本生活需求。

图5 库普工区点位图

3）反演成果。如图6，三种处理方法均能反应测线地下电阻率分布情况：采用美国某大学的REBOCC反演程序侧重于地下介质的均匀性，对于地质体整体的分布特征把握较好，处理程度较圆滑，但对于细节反应不明显；而采用经过本程序格式转换后的SCS2D反演成果能够较好的描述整体和局部特征，特别是二维反演成果，非常好的诠释了地下地质体、断裂构造的产状、埋深以及地下水位的分布特征。

图6 不同处理方法反演成果图

4）总结。通过本次工作，查清了已知水井缺水的原因，提出了更利于储水的307点；同时也推断由于深部下游水源通道的存在，下游不会因为已知水井断水而没有开发前景，为本区水源开发利用奠定了较好的基础。

4 结束语

Matlab工作平台的使用不仅大大优化了不同数据格式的转换流程，提高了野外的工作效率，也避免了优秀反演软件无法通用的弊端，多用于信号和图像处理、通讯、控制系统设计、测试和测量等众多应用领域，对地学数据处理解释起到了推动作用。

参考文献

[1]邓明，谭捍东，胡建德，等.超宽频带大地电磁信号的采集方法和技术[J].现代地质，1997，11（3）.

[2]魏文博.我国大地电磁测深新进展及瞻望[J].地球物理学进展，2002，17（2）.

[3]TAO Shi long，WAN Tian feng.Introduction of Earth Science[M].Beijing： Geological Publishing House， 1999.

[4]http：//baike.baidu.com/view/10598.htm.

作者简介

李美英（1981-），河北涿州人，助理工程师，主要从事仪器设备维修、维护及数据处理等工作。endprint

图3 运行程序调用参数文件

3）软件流程。野外EH4仪器数据采集完毕，测量台班会记录并形成一个关于点号、坐标、高程、文件号的索引参数文件，根据该文件读取原始z文件中xy方向的相关系数、视电阻率、相位、频率等，Matlab程序会根据提供的各种数据组合成SCS2D反演软件所需的mtd格式文件，并最终输出。软件大致流程如图4。

图4 软件流程图

4）程序代码及注释。

clear all；clc；

[cfilename，cdir]=uigetfile（{'*.par'，'parameter file（*.par）'；'*.*'，'All files（*.*）'}，...

'Choose a parameter file to read'）；%选择参数文件并读取。

cfile=[cdir，cfilename]；

fid=fopen（cfile，'r'）；

data=fscanf（fid，'%f'，[5，inf]）；%参数文件分四列，为点号，坐标x，坐标y，高程，文件号。

data=data'；

nsite=size（data，1）；

j=1；

for isite=1：nsite%进入循环语句。

sname=num2str（data（isite，5））；

slength=length（sname）；

if slength==1

sname=['00'，sname]；

elseif slength==2

sname=['0'，sname]；%确定z文件名。

end

zfile=['zhs.'，sname]；

zfid=fopen（zfile，'r'）；

zdata=fscanf（zfid，'%f，[15，inf]）；%读取z文件

zdata=zdata'；

for i=1：39；

if zdata（i，3）<=0；

continue；

else

mdata（j，6）=zdata（i，3）；

mdata（j，5）=zdata（i，1）；

mdata（j，1：4）=data（isite，1：4）；

mdata（j，7）=0.5；

mdata（j，8）=zdata（i，4）*3.1415926*1000/180；

mdata（j，9）=0.5；

mdata（j，10：13）=0；

mdata（j，14）=mdata（j，6）；

mdata（j，15）=mdata（j，8）；

mdata（j，16：17）=0；%生成mtd文件矩阵

j=j+1；

end

end

fclose（zfid）；

end

fclose （fid）；

mtdfile=['line'，num2str（data（1，2）），'.mtd']；

fid=fopen（mtdfile，'w'）；

fprintf（fid，'%s＼n'，'"Stn"，"GridE"，"GridN"，"Elev"，"Freq"，"ARTMobs"，"ARTMerr"，"ZPTMobs"，"ZPTMerr"，"ARTEobs"，"ARTEerr"，"ZPTEobs"，"ZPTEerr"'）；

fprintf（fid，'%s＼n'，'"From S2D2D v2.10l Date：16/12/09 Time：13：34：14"'）；

for i=1：j-1

fprintf（fid，'%6.1f，%6.1f，%6.1f，%6.1f，%6.1f，%6.1f，%6.1f，%6.1f，%6.1f'，mdata（i，1：9））；

fprintf（fid，'%s'，'，，，，'）；

fprintf（fid，'%6.1f，%6.1f'，mdata（i，14：15））；

fprintf（fid，'%s＼n'，'，，，'）；%输出mtd文件

end

fclose（fid）；

3 实例分析

1）工区概况。库普边防派出所位于奇台县八一牧场西北侧，距奇台县城约197千米，为昌吉边防支队的下属单位，多年来一直处于缺水状态。为做好双拥工作，更为解决官兵的实际用水困难，应昌吉国土资源局邀请，物探队派遣EH4项目组，了解该区地下电阻率分布特征，为寻找储水构造并布置钻井提供依据。

该区地处东准噶尔碱性花岗岩体中，寻找裂隙水是本次工作的主攻目标。

2）工作布置。如图5，限于工区实际条件，布置测线：3线、4线、5线。取点数较多的3线为例，3线307点西北侧为学校办公楼，楼前308点西约6米处为已知水井，该水井为前些年所打，井深10米，基底为花岗岩。由于上游水量降低，该井出水量日益减少，已不能满足学校最基本生活需求。

图5 库普工区点位图

3）反演成果。如图6，三种处理方法均能反应测线地下电阻率分布情况：采用美国某大学的REBOCC反演程序侧重于地下介质的均匀性，对于地质体整体的分布特征把握较好，处理程度较圆滑，但对于细节反应不明显；而采用经过本程序格式转换后的SCS2D反演成果能够较好的描述整体和局部特征，特别是二维反演成果，非常好的诠释了地下地质体、断裂构造的产状、埋深以及地下水位的分布特征。

图6 不同处理方法反演成果图

4）总结。通过本次工作，查清了已知水井缺水的原因，提出了更利于储水的307点；同时也推断由于深部下游水源通道的存在，下游不会因为已知水井断水而没有开发前景，为本区水源开发利用奠定了较好的基础。

4 结束语

Matlab工作平台的使用不仅大大优化了不同数据格式的转换流程，提高了野外的工作效率，也避免了优秀反演软件无法通用的弊端，多用于信号和图像处理、通讯、控制系统设计、测试和测量等众多应用领域，对地学数据处理解释起到了推动作用。

参考文献

[1]邓明，谭捍东，胡建德，等.超宽频带大地电磁信号的采集方法和技术[J].现代地质，1997，11（3）.

[2]魏文博.我国大地电磁测深新进展及瞻望[J].地球物理学进展，2002，17（2）.

[3]TAO Shi long，WAN Tian feng.Introduction of Earth Science[M].Beijing： Geological Publishing House， 1999.

[4]http：//baike.baidu.com/view/10598.htm.

作者简介

李美英（1981-），河北涿州人，助理工程师，主要从事仪器设备维修、维护及数据处理等工作。endprint

图3 运行程序调用参数文件

3）软件流程。野外EH4仪器数据采集完毕，测量台班会记录并形成一个关于点号、坐标、高程、文件号的索引参数文件，根据该文件读取原始z文件中xy方向的相关系数、视电阻率、相位、频率等，Matlab程序会根据提供的各种数据组合成SCS2D反演软件所需的mtd格式文件，并最终输出。软件大致流程如图4。

图4 软件流程图

4）程序代码及注释。

clear all；clc；

[cfilename，cdir]=uigetfile（{'*.par'，'parameter file（*.par）'；'*.*'，'All files（*.*）'}，...

'Choose a parameter file to read'）；%选择参数文件并读取。

cfile=[cdir，cfilename]；

fid=fopen（cfile，'r'）；

data=fscanf（fid，'%f'，[5，inf]）；%参数文件分四列，为点号，坐标x，坐标y，高程，文件号。

data=data'；

nsite=size（data，1）；

j=1；

for isite=1：nsite%进入循环语句。

sname=num2str（data（isite，5））；

slength=length（sname）；

if slength==1

sname=['00'，sname]；

elseif slength==2

sname=['0'，sname]；%确定z文件名。

end

zfile=['zhs.'，sname]；

zfid=fopen（zfile，'r'）；

zdata=fscanf（zfid，'%f，[15，inf]）；%读取z文件

zdata=zdata'；

for i=1：39；

if zdata（i，3）<=0；

continue；

else

mdata（j，6）=zdata（i，3）；

mdata（j，5）=zdata（i，1）；

mdata（j，1：4）=data（isite，1：4）；

mdata（j，7）=0.5；

mdata（j，8）=zdata（i，4）*3.1415926*1000/180；

mdata（j，9）=0.5；

mdata（j，10：13）=0；

mdata（j，14）=mdata（j，6）；

mdata（j，15）=mdata（j，8）；

mdata（j，16：17）=0；%生成mtd文件矩阵

j=j+1；

end

end

fclose（zfid）；

end

fclose （fid）；

mtdfile=['line'，num2str（data（1，2）），'.mtd']；

fid=fopen（mtdfile，'w'）；

fprintf（fid，'%s＼n'，'"Stn"，"GridE"，"GridN"，"Elev"，"Freq"，"ARTMobs"，"ARTMerr"，"ZPTMobs"，"ZPTMerr"，"ARTEobs"，"ARTEerr"，"ZPTEobs"，"ZPTEerr"'）；

fprintf（fid，'%s＼n'，'"From S2D2D v2.10l Date：16/12/09 Time：13：34：14"'）；

for i=1：j-1

fprintf（fid，'%6.1f，%6.1f，%6.1f，%6.1f，%6.1f，%6.1f，%6.1f，%6.1f，%6.1f'，mdata（i，1：9））；

fprintf（fid，'%s'，'，，，，'）；

fprintf（fid，'%6.1f，%6.1f'，mdata（i，14：15））；

fprintf（fid，'%s＼n'，'，，，'）；%输出mtd文件

end

fclose（fid）；

3 实例分析

1）工区概况。库普边防派出所位于奇台县八一牧场西北侧，距奇台县城约197千米，为昌吉边防支队的下属单位，多年来一直处于缺水状态。为做好双拥工作，更为解决官兵的实际用水困难，应昌吉国土资源局邀请，物探队派遣EH4项目组，了解该区地下电阻率分布特征，为寻找储水构造并布置钻井提供依据。

该区地处东准噶尔碱性花岗岩体中，寻找裂隙水是本次工作的主攻目标。

2）工作布置。如图5，限于工区实际条件，布置测线：3线、4线、5线。取点数较多的3线为例，3线307点西北侧为学校办公楼，楼前308点西约6米处为已知水井，该水井为前些年所打，井深10米，基底为花岗岩。由于上游水量降低，该井出水量日益减少，已不能满足学校最基本生活需求。

图5 库普工区点位图

3）反演成果。如图6，三种处理方法均能反应测线地下电阻率分布情况：采用美国某大学的REBOCC反演程序侧重于地下介质的均匀性，对于地质体整体的分布特征把握较好，处理程度较圆滑，但对于细节反应不明显；而采用经过本程序格式转换后的SCS2D反演成果能够较好的描述整体和局部特征，特别是二维反演成果，非常好的诠释了地下地质体、断裂构造的产状、埋深以及地下水位的分布特征。

图6 不同处理方法反演成果图

4）总结。通过本次工作，查清了已知水井缺水的原因，提出了更利于储水的307点；同时也推断由于深部下游水源通道的存在，下游不会因为已知水井断水而没有开发前景，为本区水源开发利用奠定了较好的基础。

4 结束语

Matlab工作平台的使用不仅大大优化了不同数据格式的转换流程，提高了野外的工作效率，也避免了优秀反演软件无法通用的弊端，多用于信号和图像处理、通讯、控制系统设计、测试和测量等众多应用领域，对地学数据处理解释起到了推动作用。

参考文献

[1]邓明，谭捍东，胡建德，等.超宽频带大地电磁信号的采集方法和技术[J].现代地质，1997，11（3）.

[2]魏文博.我国大地电磁测深新进展及瞻望[J].地球物理学进展，2002，17（2）.

[3]TAO Shi long，WAN Tian feng.Introduction of Earth Science[M].Beijing： Geological Publishing House， 1999.

[4]http：//baike.baidu.com/view/10598.htm.

作者简介

李美英（1981-），河北涿州人，助理工程师，主要从事仪器设备维修、维护及数据处理等工作。endprint