EH4数据处理在SCS2D二维反演平台中的实现
作者:李美英 屈栓柱
来源:《硅谷》2014年第05期
        摘 要 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能够深入到科学研究及工程计算各个领域的重要原因。
        图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)
        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));

版权声明:本站内容均来自互联网,仅供演示用,请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198,我们将在24小时内删除。