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