一種基于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-），男，本科，助理工程師，研究方向：物探數據處理。