Matlab在联剖法数据处理中的应用

聂小力, 许瑞华, 李宇杰, 罗敏玄, 谌旭辉

(1.中国地质大学 地球物理与信息技术学院，北京 100083；2.武警黄金第十一支队，拉萨 850000)

Matlab在联剖法数据处理中的应用

聂小力1，2, 许瑞华2, 李宇杰2, 罗敏玄2, 谌旭辉2

(1.中国地质大学 地球物理与信息技术学院，北京 100083；2.武警黄金第十一支队，拉萨 850000)

考虑到Matlab编程效率高、使用方便等优点，以联合剖面法数据处理为例，分别利用Matlab，编写了联剖法数据处理过程中的地形校正、单线剖面图绘制、剖平图绘制、点位投影等程序,实现了对联剖法数据的全套处理，利用西藏某矿区实测数据给出了处理实例。利用Matlab图形用户界面(GUI)设计编写了可视化的联剖法数据处理程序，实现了联剖法数据处理的可视化操作，简单、便捷。

MATLAB； 联合剖面法； 地形校正； 图形用户界面GUI

0 前言

地球物理数据处理是物探工作中的重要步骤，而直流电法的数据处理目前尚没有较为全面的成型软件，数据预处理、地形改正、各电性参数的计算以及成图往往需要调用不同的软件，软件质量及处理效果也参差不齐。这使得数据处理的过程系统性不强，较为盲目。而Matlab作为在国际上被广泛接受和使用的数值计算软件，将其运用到物探数据处理当中，每个处理环节的工作都可以很好的完成，使得整个处理过程系统性强、对比性高，让处理后期的异常解释更为直观可靠。

这里利用Matlab,以联剖法数据处理为例，通过编写程序基本完成了联剖法数据的全套处理，并且利用图形用户界面GUI编写了一套可视化的联剖法数据处理软件，使处理过程更为简洁直观[1]。也希望以此为窗口，启发更多地质工作者将Matlab运用到地质领域的各个方面中去。

1 联剖法中的地形校正

1.1 地形校正计算方法

在电法勘探中,起伏的地形往往使有用异常面目皆非。因而在对电法数据资料进行地质解释时, 克服地形影响一直是重要问题之一。为克服地形影响，需要计算求出由起伏地形引起的纯地形影响 。这里利用Matlab,通过数值模拟来计算纯地形影响。

选取的地改方法较为简单。即利用地表与水平面之间的夹角α，以及供电电极与记录点的高程数据，根据视电阻率的微分表达式(式(1))[2]，计算求取纯地形影响。

(1)

jMN*H′=jO*H

(2)

由此可得

(3)

进而求得归一化的纯地形异常式(4)。

(4)

利用纯地形异常即可根据式(4)求得进过地改后的视电阻率值：

(5)

式中：ρS为野外实测电阻率；AO′、sin∠O′AO可根据供电电极、记录点的位置及高程求得。

图1 地形改正计算模型Fig.1 Calculation model for topographic correction

1.2 利用Matlab对不同地形进行模拟计算的实例

图2是用纯地形影响求解方法求取的二维山脊和二维山谷地形上的联合剖面视电阻率曲线。

图2 模拟计算模型及纯地形影响剖面图Fig.2 Relief model and the result of pure topographic influence(a)山脊地形上联剖剖面曲线； (b)山谷地形上联剖剖面曲线

2 单线剖面图绘制

在分析单条地电勘探线数据时，需要先绘制单线剖面图。绘制剖面图时可选取不同电性参数，从而从不同角度对地电测量结果进行分析。这其中能反映测线电性特征的参数主要有以下几种：①视电阻率；②视极化率；③比值参数；④岐离带参数等。

在作者编写的单线剖面图Matlab程序[4]中，可以选择是否进行地形校正，并可以选取需要输出的参数类型或者输出全部参数曲线。下面根据实测数据给出具体实例。

2.1 地形校正

在做剖面图时，当地形起伏比较大时，会对测量结果产生较大影响，所以在绘制图形前需要确定是否需要进行地改处理[5]。

HIGH=xlsread('data.xlsx'); %读取工区高程数据

circle=(AO-step/2)/step; %求AO间相隔的点数

high=HIGH((n-circle):(k+circle),line+2);

for i=1:1:size(m); %求地形起伏角度

ii=i+circle;

angleA(i)=(high(i)-high(ii))/step/circle;

angleB(i)=(high(ii)-high(ii+circle))/step/circle;

远程脉冲水表的动力装置是叶轮，当水流通过叶轮时，带动叶轮旋转，而水流的流速与叶轮的转速成正比，因水流驱动叶轮出喷口的截面积为常数，故叶轮的转速与流量也成正比。因生产用水中含有少量杂质，长时间运行后，水表叶轮会逐步聚集缠绕部分纤维物，导致水表不能正常工作。

end

for j=1:1:size(angleA,2); %求取纯地形影响

corValueA(j)=1/(1-angleA(j));

corValueB(j)=1/(1+angleB(j));

end

for kk=1:1:size(angleA,2); %求地改后视电阻率

psAC(kk)=dataA(kk)/corValueA(kk);

psBC(kk)=dataB(kk)/corValueB(kk);

end

图3中，观察地形可以发现，水平距离200 m的高差达到了50 m，这将对测量结果造成影响。

对比图3(a),图3(b)可以发现，曲线发生了变化。①从图3(a)中可以看到，视电阻率曲线245 m处存在一个低阻正交点，而从图3(b)中可以看到，这一交点其实是不存在的，而是受地形起伏影响而产生的“假交点”。②图3(a)中，410 m位置处存在一个低阻正交点，在图3(b)中，这一交点出现在370 m处，这说明地形起伏的影响使得交点的位置发生了变化。

图3 地形校正曲线图Fig.3 Profile of topographic correction(a)野外实测视电阻率剖面图；(b)地改后的视电阻率剖面图；(c)纯地形影响剖面图；(d)地形剖面图

由图3可以看出，当地形起伏较大时，倘若不对原始数据进行地形校正，直接根据交点性质与位置加以解释的话，会得到错误的推断。

2.2 比值参数

当观测结果受地表局部不均匀体影响时，观测到的联合剖面曲线将出现同步跳跃现象。这将会给基于交点性质来解释观测结果的联合剖面法造成一定的困难。为此，可以通过计算比值参数来对地表不均匀体的影响进行压制。比值参数按式(6)进行计算。

(6)

求取比值参数和岐离带参数的程序代码如下：

Fa=psAC./psBC; %计算求取Fa

Fb=psBC./psAC; %计算求取Fb

subplot(2,1,1); %分幅绘制比值曲线

plot(x,Fa,'color','r');

hold on;

plot(x,Fb,'color','b');

for i=1:1(size(psAC)-1); %计算岐离带参数

lamutaFa=Fa(i)/Fa(i+1);

end

subplot(2,1,2); %分幅绘制岐离带参数曲线

plot(x,lamuta,'r');

观察图4可以看出，与经过地形改正后的视电阻率曲线进行比较，虽然地形改正消除了曲线受地形的影响，但由于局部不均匀体的影响，视电阻率曲线靠的更近，同时同步跳跃。计算比值参数后，比值曲线分离较大，异常的交点更加明显。

2.3 歧离带参数

歧离带参数是一种用来突出地质异常的比值参数。对于低阻直立薄板上，联合剖面的歧离带参数曲线为单峰异常。模型模拟结果表明，该参数具有一定的复合异常分析能力。式(7)为岐离带参数的计算公式。

(7)

图5是利用野外实测数据通过计算得到的岐离带参数曲线。在视电阻率曲线中，低阻正交点往往指示着低阻体的存在，而对比观察岐离带曲线与视电阻率曲线可以看出，剖面290 m、435 m、700 m、920 m出现的低阻正交点都有相应的岐离带单峰异常出现，这说明岐离带参数对于低阻体的存在确实有着指示作用。但需要注意的是，这一参数受干扰影响大，在实际运用中需谨慎使用。

图4 比值参数曲线图Fig.4 Profile of ratio parameters(a)视电阻率曲线；(b)比值参数曲线

图5 岐离带参数曲线图Fig.5 Profile of straggling parameters(a)视电阻率曲线；(b)岐离带参数曲线

3 平面剖面图绘制

平面剖面图是直流电法勘探中一种可以同时显示某种地球物理场剖面及其平面特征的重要图件。

平剖图中，各测线异常都用剖面形式表示，且各剖面比例尺统一。因而通过各测线之间的对比分析能清楚地表示异常平面和剖面的分布规律。进一步探明矿区异常的走向及其延展性。

在作者编写的Matlab平剖图绘制程序中，可自行选择需要输出的剖面号，选择是否进行地形校正，并且可以选择需要输出的参数类型(包括视电阻率、视极化率、比值参数、岐离带参数)。程序输出图形，见图6。对于这些输出参数的控制，具体的Matab程序代码如下：

%以矩阵形式输出所需处理的测线号

%选择是否地改并选择需要输出的曲线类型

inputlines=input('请输入处理线号：');

cor=input('是否地改(YES:1,NO:2)');

CurveType=input('输出曲线类型(视电阻率：1 比值曲线：2 岐离带曲线：3)')

%根据输入的线号编程循环输出各剖面曲线

图6 不同参数的平剖图曲线Fig.6 Profile of different parameters(a)带地改的视电阻率曲线；(b)不带地改的比值参数曲线

4 GUI设计的可视化程序

为使数据处理及成图过程更为直观、便捷，利用Matlab GUI设计了以下可视化程序。该程序共有三个处理模块：①单线剖面图绘制；②平剖图绘制；③点位投影。图7为编写的主菜单界面，用以选择调用处理模块。

图7 联剖法数据处理程序主菜单Fig.7 Main menu of the data processing software for composite profiling method

单线剖面绘制程序用以输出单条地电断面的剖面图。程序首先需要选择输入Excel格式的数据文件及输出的线号，即可输出剖面图。做图过程中可以随时更改输出线号，输出参数，是否地改，输出剖面的范围等。此外还可以设置图表输出参数来使图像美观。程序界面见图8。

平剖图绘制程序用以绘制平面剖面图。选择Excel格式的数据文件，并输入需要输出的所有测线号，即可绘制出平剖图。同时可以选择输出的电性参数，是否地改，图像会随设定参数的改变随时更新。其次还能设置图表参数。程序界面见图9。

点位投影用来输出测区的三维地形，并利用视电阻率联剖数据计算求得对称四极装置的视电阻率异常，将其投影到三维地形上并上色。其次还将投影所有测点位置，以及程序拾取得到的正反交点位置到三维地形上，以便更为直观的分析观察测量结果。程序界面见图10。

图8 单线剖面图绘制程序Fig.8 Mapping program for single profile

图9 平剖图绘制程序Fig.9 Mapping program for profile map

5 结束语

通过Matlab编程基本实现了联剖法数据整理及成图的全套处理。而GUI编写的可视化程序使得处理过程简单直观。当然所编写的程序中地形改正方法还有待进一步改善，数据文件的兼容性问题也需要进一步提高。总之Matlab强大的数据处理功能使得它在整个地学领域都有其用武之处。希望能看到Matlab在地质工作中得到更为广泛地使用。

[1] 张志涌. 精通Matlab R2011a [M]. 北京：北京航天航空大学出版社，2011. ZHANG Z Y. Proficient in Matlab R2011a [M]. Beijing：Beihang University Press, 2011. (In Chinese)

[2] 傅良魁. 应用地球物理教程—电法 放射性 地热 [M]. 北京：地质出版社，1991. FU L K. Applied geophysics [M]. Beijing : Geological Publishing House, 1991. (In Chinese)

图10 点位投影程序Fig.10 Mapping program for point pojection

[3] 李金铭. 地电场与电法勘探 [M]. 北京：地质出版社, 2005． LI J M. Electric and electrical prospecting [M].Beijing:Geological Publishing House,2005. (In Chinese)

[4] 童孝忠,王涛,柳建新. Matlab软件包在地球物理数值计算与辅助教学中的应用[J].物探化探计算技术, 2013,35(5): 512-518. TONG X Z, WANG T, LIU J X. Numerical calculation and aided teaching in geophysics using Matlab [J]. Computing Techniques for Geophysics and Geochemical exploration, 2013,35(5):512-518. (In Chinese)

[5] 汤井田,辛会翠,王冉. 点电源下复杂角域地形影响及校正[J]. 吉林大学学报(地球科学版),2012, 42(1): 254—261. TANG J T, XIN H C, WANG R. Terrain effect of complicated angle-damain and topographic correction for a point electric source [J].Journal of Jilin University( Earth Science Edition), 2012, 42(1):254-261. (In Chinese)

[6] 许瑞华. 比值参数在激电联剖视极化率数据处理中的应用[J]. 四川地质学报,2012, 32(增刊): 72-74. XU R H. Data processing for composite profiling method using ratio parameters[J]. Acta Geologica Sichuan, 2012, 32(S):72-74. (In Chinese)

[7] 陈义群, 陈华. 基于MATLAB的工程物探软件快速开发[J]. 地球物理学进展,2004, 19(4): 802-806. CHEN Y Q, CHEN H. Fast developmet of engineering geophysicl software based on Matlab[J]. Progress In Geophysics, 2004, 19(4):802-806. (In Chinese)

[8] 王鹏.基于MATLAB的煤矿TEM数据体三维可视化技术[J].地球物理学进展,2014,29(3):1277-1283. WANG P.Three-dimensional visualization of colliery TEM data base on Matlab[J].Progress In Geophysics, 2014, 29(3):1277-1283. (In Chinese)

Data processing for composite profiling method using Matlab

NIE Xiaoli1,2， XU Ruihua2， LI Yujie2， LUO Minxuan2， CHENG Xuhui2

(1. Hina University of Geoscience,School of Imfomation Technology and Geophysics,Beijing 100083，China；2. No.11 Gold Geological Party of the CAPF, Lhasa 850000,China)

Matlab programming is easy to use and with high efficiency. Considering its advantages, using MATLAB in data processing of composite profiling method as a case, and writing programs for single profile mapping, profile map drawing, point projecting, basically complete the full set of data processing of the method. Secondly, using MATLAB graphical user interface (GUI) design, writing the visual data processing program of composite profiling method.

Matlab; composite profiling method; topographic correction; graphical user interface(GUI)

2015-11-17 改回日期：2015-12-24

中国地质调查项目(12120114083501)

聂小力(1991-)，男，助理工程师，从事地球物理勘探方法研究， E-mail：niexiaoli0912@163.com。

1001-1749(2017)01-0137-07

P 631.4

A

10.3969/j.issn.1001-1749.2017.01.20