楊凱

（中國人民武裝警察部隊黃金第二支隊，內(nèi)蒙古呼和浩特 010010）

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

楊凱

（中國人民武裝警察部隊黃金第二支隊，內(nèi)蒙古呼和浩特 010010）

測深是勘察地下電阻率或其他物理特征分布的重要地球物理方法，主要分為電阻率測深、激電測深、電磁測深等，一般只需得到測深剖面的一二維斷面圖就可達(dá)到一定的勘探目的，但有些時候為了從不同角度、不同側(cè)重點來研究異常，還要有不同深度的水平電阻率等值線圖?；诖耍ㄟ^FORTRAN語言編程，利用二維插值原理來統(tǒng)一提取多條測深剖面相應(yīng)深度的電阻率數(shù)據(jù)，進(jìn)而繪制水平電阻率等值線圖。

Fortran 90；測深；二維插值；水平切面

FORTRAN語言是世界上最早出現(xiàn)的高級程序設(shè)計語言，其以嚴(yán)謹(jǐn)、規(guī)范、高效的特點普遍應(yīng)用于數(shù)值分析、工程計算和科學(xué)研究等領(lǐng)域［1］。利用FORTRAN可以對物探數(shù)據(jù)進(jìn)行批量、高效的處理。本文就是筆者在實踐工作中探索出的根據(jù)二維差值原理利用FORTRAN編程實現(xiàn)從多條測深剖面數(shù)據(jù)中提取相應(yīng)深度水平切面的過程。

1 基本原理

根據(jù)給定矩形域上n×m個結(jié)點上的函數(shù)值，利用二元插值公式計算出指定插值點處的函數(shù)值。

已知矩形域上n×m個結(jié)點的2個方向上的坐標(biāo)分別為：

其相應(yīng)的函數(shù)值為：

計算插值點（u，v）為中心，在X方向上，前后各取4個坐標(biāo)：

在Y方向上，前后也取4個坐標(biāo)：

然后利用二元插值公式：

計算插值點（u，v）處的函數(shù)值。

根據(jù)以上原理即可通過二元插值在電阻率等值斷面圖GRD網(wǎng)格化文件中提取出所需深度的電阻率值。

2 程序流程

程序處理流程主要包括測深剖面GRD網(wǎng)格化數(shù)據(jù)輸入，用戶輸入相應(yīng)參數(shù)，計算電阻率值并輸出，以下為主要模塊的程序代碼。

2.1 數(shù)據(jù)輸入模塊

該功能以模塊化子例行程序?qū)崿F(xiàn)，命名為input_da?ta2，其中inputname1為輸入GRD網(wǎng)格化文件名，Z是一個二維數(shù)組用于存放網(wǎng)格化文件數(shù)據(jù)，N和M分別為橫向和縱向的數(shù)據(jù)個數(shù)，xmin1和xmax1為橫向起始坐標(biāo)值和結(jié)束坐標(biāo)值，ymin1和ymax1為縱向起始坐標(biāo)值和結(jié)束坐標(biāo)值，以下為程序代碼：

2.2 參數(shù)輸入模塊

處理參數(shù)主要從兩部分輸入，一是CMD輸入輸出控制文件，其功能主要是存放待處理的測深剖面GRD網(wǎng)格化文件名，二是主程序中的輸入?yún)?shù)設(shè)置，主要方式以人機(jī)對話為主，主要輸入測線數(shù)L、測線方位（設(shè)定x為東西向，y為南北向）和線距LD以及截取深度標(biāo)高h(yuǎn)eight、數(shù)據(jù)輸出文件名outputname等，代碼如下：

2.3 電阻率值提取模塊

該功能以模塊化子例行程序?qū)崿F(xiàn)，命名為interpola?tion，其借鑒了徐世良［2］先生FORTRAN常用算法集中的二維差值算法ESLGQ（X，Y，Z，N，M，U，V，W）子例行程序，其中WW為存放提取電阻率值的一維數(shù)組，XX與YY為存放水平切面橫縱向坐標(biāo)的一維數(shù)組，height實形變量表征為所要截取的深度標(biāo)高，該模塊只實現(xiàn)了單個測線的處理，要實現(xiàn)多條測線處理還得循環(huán)調(diào)用該模塊，實現(xiàn)測線的批量處理，程序代碼如下：

process為計算GRD網(wǎng)格化文件中的數(shù)據(jù)節(jié)點坐標(biāo)值的子例行程序，該坐標(biāo)值用于差值計算。代碼如下：

2.4 計算輸出模塊

該段程序內(nèi)包含有2個子例行程序input_NM和in?terpolation，interpolation已在上文中介紹，input_NM的主要功能是從每條測線GRD網(wǎng)格化文件中提取N和M即橫向和縱向的數(shù)據(jù)個數(shù)，通過循環(huán)調(diào)用這2個子例行程序，就可以得到每條測線的電阻率值，并以x，y，z三列數(shù)據(jù)的格式寫入到用戶指定文件中，再經(jīng)過網(wǎng)格化或其他處理手段可以得到相應(yīng)的圖件和成果，為找礦提供依據(jù)。程序代碼如下：

通過編寫主程序代碼，分別按照處理步驟依次調(diào)用各個功能模塊即可實現(xiàn)數(shù)據(jù)的處理，本文中的代碼均在FORTRAN90中檢驗通過。

3 實際應(yīng)用

以內(nèi)蒙古自治區(qū)東烏旗某礦區(qū)所做的測深工作為應(yīng)用對象，測區(qū)共開展可控源音頻大地電磁測深剖面11條，測線長1 640m，點距40m，線距200m。測區(qū)的高阻體主要為凝灰?guī)r、花崗巖等巖體，地層主要以中低阻為主要特征，控礦構(gòu)造主要以低電阻率特征的斷裂帶和接觸帶為主，但該地區(qū)的地表覆蓋較厚，一些隱伏的控礦構(gòu)造利用常規(guī)方法不易揭露，為了更好地掌握深部電阻率的分布規(guī)律，利用該方法對測區(qū)的11條測深剖面進(jìn)行了800、850、900、950m 4個深度標(biāo)高的電阻率提取成圖處理（見圖1），通過該處理結(jié)果很好地掌握了測區(qū)不同深度上的電阻率分布特征，圈定了成礦有利部位，為下一步工作開展提供了重要依據(jù)［3］。

圖1 測區(qū)深度標(biāo)高800、850、900、950m水平面聯(lián)合電阻率等值線圖

4 結(jié)論

①常規(guī)的二維差值方法只能在一條測深剖面等值線圖上提取相應(yīng)深度的數(shù)據(jù)，為了繪制水平等值線圖必須有多個測深剖面的數(shù)據(jù)，筆者利用FORTRAN處理數(shù)據(jù)的優(yōu)點通過子例行程序調(diào)用和循環(huán)結(jié)構(gòu)解決了同時提取多條剖面中同一深度電阻率數(shù)據(jù)的問題，簡化了工作步驟，提高了工作效率。

②通過實際應(yīng)用可以看到，該方法所成圖件清晰準(zhǔn)確，很好地反映了異常在平面上的分布特征，為從不同角度研究異常提供了條件，同時也檢驗了該方法的正確性和穩(wěn)定性。

［1］白云.FORTRAN90程序設(shè)計［M］.上海：華東理工大學(xué)出版社，2003.

［2］徐世良.FORTRAN常用算法集［M］.北京：清華大學(xué)出版社，1995.

［3］何繼善.可控源音頻大地電磁法［M］.長沙：中南工業(yè)大學(xué)出版社，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-），男，本科，助理工程師，研究方向：物探數(shù)據(jù)處理。