一种基于FORTRAN的多测深剖面水平切面提取方法

杨凯

（中国人民武装警察部队黄金第二支队，内蒙古呼和浩特 010010）

一种基于FORTRAN的多测深剖面水平切面提取方法

杨凯

（中国人民武装警察部队黄金第二支队，内蒙古呼和浩特 010010）

测深是勘察地下电阻率或其他物理特征分布的重要地球物理方法，主要分为电阻率测深、激电测深、电磁测深等，一般只需得到测深剖面的一二维断面图就可达到一定的勘探目的，但有些时候为了从不同角度、不同侧重点来研究异常，还要有不同深度的水平电阻率等值线图。基于此，通过FORTRAN语言编程，利用二维插值原理来统一提取多条测深剖面相应深度的电阻率数据，进而绘制水平电阻率等值线图。

Fortran 90；测深；二维插值；水平切面

FORTRAN语言是世界上最早出现的高级程序设计语言，其以严谨、规范、高效的特点普遍应用于数值分析、工程计算和科学研究等领域［1］。利用FORTRAN可以对物探数据进行批量、高效的处理。本文就是笔者在实践工作中探索出的根据二维差值原理利用FORTRAN编程实现从多条测深剖面数据中提取相应深度水平切面的过程。

1 基本原理

根据给定矩形域上n×m个结点上的函数值，利用二元插值公式计算出指定插值点处的函数值。

已知矩形域上n×m个结点的2个方向上的坐标分别为：

其相应的函数值为：

计算插值点（u，v）为中心，在X方向上，前后各取4个坐标：

在Y方向上，前后也取4个坐标：

然后利用二元插值公式：

计算插值点（u，v）处的函数值。

根据以上原理即可通过二元插值在电阻率等值断面图GRD网格化文件中提取出所需深度的电阻率值。

2 程序流程

程序处理流程主要包括测深剖面GRD网格化数据输入，用户输入相应参数，计算电阻率值并输出，以下为主要模块的程序代码。

2.1 数据输入模块

该功能以模块化子例行程序实现，命名为input_da⁃ta2，其中inputname1为输入GRD网格化文件名，Z是一个二维数组用于存放网格化文件数据，N和M分别为横向和纵向的数据个数，xmin1和xmax1为横向起始坐标值和结束坐标值，ymin1和ymax1为纵向起始坐标值和结束坐标值，以下为程序代码：

2.2 参数输入模块

处理参数主要从两部分输入，一是CMD输入输出控制文件，其功能主要是存放待处理的测深剖面GRD网格化文件名，二是主程序中的输入参数设置，主要方式以人机对话为主，主要输入测线数L、测线方位（设定x为东西向，y为南北向）和线距LD以及截取深度标高height、数据输出文件名outputname等，代码如下：

2.3 电阻率值提取模块

该功能以模块化子例行程序实现，命名为interpola⁃tion，其借鉴了徐世良［2］先生FORTRAN常用算法集中的二维差值算法ESLGQ（X，Y，Z，N，M，U，V，W）子例行程序，其中WW为存放提取电阻率值的一维数组，XX与YY为存放水平切面横纵向坐标的一维数组，height实形变量表征为所要截取的深度标高，该模块只实现了单个测线的处理，要实现多条测线处理还得循环调用该模块，实现测线的批量处理，程序代码如下：

process为计算GRD网格化文件中的数据节点坐标值的子例行程序，该坐标值用于差值计算。代码如下：

2.4 计算输出模块

该段程序内包含有2个子例行程序input_NM和in⁃terpolation，interpolation已在上文中介绍，input_NM的主要功能是从每条测线GRD网格化文件中提取N和M即横向和纵向的数据个数，通过循环调用这2个子例行程序，就可以得到每条测线的电阻率值，并以x，y，z三列数据的格式写入到用户指定文件中，再经过网格化或其他处理手段可以得到相应的图件和成果，为找矿提供依据。程序代码如下：

通过编写主程序代码，分别按照处理步骤依次调用各个功能模块即可实现数据的处理，本文中的代码均在FORTRAN90中检验通过。

3 实际应用

以内蒙古自治区东乌旗某矿区所做的测深工作为应用对象，测区共开展可控源音频大地电磁测深剖面11条，测线长1 640m，点距40m，线距200m。测区的高阻体主要为凝灰岩、花岗岩等岩体，地层主要以中低阻为主要特征，控矿构造主要以低电阻率特征的断裂带和接触带为主，但该地区的地表覆盖较厚，一些隐伏的控矿构造利用常规方法不易揭露，为了更好地掌握深部电阻率的分布规律，利用该方法对测区的11条测深剖面进行了800、850、900、950m 4个深度标高的电阻率提取成图处理（见图1），通过该处理结果很好地掌握了测区不同深度上的电阻率分布特征，圈定了成矿有利部位，为下一步工作开展提供了重要依据［3］。

图1 测区深度标高800、850、900、950m水平面联合电阻率等值线图

4 结论

①常规的二维差值方法只能在一条测深剖面等值线图上提取相应深度的数据，为了绘制水平等值线图必须有多个测深剖面的数据，笔者利用FORTRAN处理数据的优点通过子例行程序调用和循环结构解决了同时提取多条剖面中同一深度电阻率数据的问题，简化了工作步骤，提高了工作效率。

②通过实际应用可以看到，该方法所成图件清晰准确，很好地反映了异常在平面上的分布特征，为从不同角度研究异常提供了条件，同时也检验了该方法的正确性和稳定性。

［1］白云.FORTRAN90程序设计［M］.上海：华东理工大学出版社，2003.

［2］徐世良.FORTRAN常用算法集［M］.北京：清华大学出版社，1995.

［3］何继善.可控源音频大地电磁法［M］.长沙：中南工业大学出版社，1990.

A Method for Extracting Horizontal Plane Section of Multi Depth Sounding Profile Based on FORTRAN

Yang Kai
（Chinese People’S Armed Police Force Gold Second Detachment，Hohhot Inner Mongolia 010010）

Sounding is an important geophysical method to survey the distribution of underground resistivity or other physical characteristics.It is mainly divided into resisitivity sounding,IP sounding,electromagnetic sounding and so on.Under normal circumstances,we can achieve prospecting purpose just to get one-dimensional or two-dimension⁃al profile contour map,but sometime we need to get resistivity horizontal layers contour map of different heights to re⁃search anomalies from different perspectives and different emphases.Based on this,through the FORTRAN language programming,the horizontal resistivity contour map was draw,by using the two-dimensional interpolation principle to extract the resistivity data of the corresponding depth of the multi depth sounding section.

fortran 90；depth sounding；two-dimensional interpolation；horizontal section

P631.3

A

1003-5168（2016）12-0058-03

2016-11-01

杨凯（1991-），男，本科，助理工程师，研究方向：物探数据处理。