Matlab在聯(lián)剖法數(shù)據(jù)處理中的應(yīng)用

聶小力, 許瑞華, 李宇杰, 羅敏玄, 諶旭輝

(1.中國(guó)地質(zhì)大學(xué) 地球物理與信息技術(shù)學(xué)院，北京 100083；2.武警黃金第十一支隊(duì)，拉薩 850000)

Matlab在聯(lián)剖法數(shù)據(jù)處理中的應(yīng)用

聶小力1，2, 許瑞華2, 李宇杰2, 羅敏玄2, 諶旭輝2

(1.中國(guó)地質(zhì)大學(xué) 地球物理與信息技術(shù)學(xué)院，北京 100083；2.武警黃金第十一支隊(duì)，拉薩 850000)

考慮到Matlab編程效率高、使用方便等優(yōu)點(diǎn)，以聯(lián)合剖面法數(shù)據(jù)處理為例，分別利用Matlab，編寫了聯(lián)剖法數(shù)據(jù)處理過(guò)程中的地形校正、單線剖面圖繪制、剖平圖繪制、點(diǎn)位投影等程序,實(shí)現(xiàn)了對(duì)聯(lián)剖法數(shù)據(jù)的全套處理，利用西藏某礦區(qū)實(shí)測(cè)數(shù)據(jù)給出了處理實(shí)例。利用Matlab圖形用戶界面(GUI)設(shè)計(jì)編寫了可視化的聯(lián)剖法數(shù)據(jù)處理程序，實(shí)現(xiàn)了聯(lián)剖法數(shù)據(jù)處理的可視化操作，簡(jiǎn)單、便捷。

MATLAB； 聯(lián)合剖面法； 地形校正； 圖形用戶界面GUI

0 前言

地球物理數(shù)據(jù)處理是物探工作中的重要步驟，而直流電法的數(shù)據(jù)處理目前尚沒(méi)有較為全面的成型軟件，數(shù)據(jù)預(yù)處理、地形改正、各電性參數(shù)的計(jì)算以及成圖往往需要調(diào)用不同的軟件，軟件質(zhì)量及處理效果也參差不齊。這使得數(shù)據(jù)處理的過(guò)程系統(tǒng)性不強(qiáng)，較為盲目。而Matlab作為在國(guó)際上被廣泛接受和使用的數(shù)值計(jì)算軟件，將其運(yùn)用到物探數(shù)據(jù)處理當(dāng)中，每個(gè)處理環(huán)節(jié)的工作都可以很好的完成，使得整個(gè)處理過(guò)程系統(tǒng)性強(qiáng)、對(duì)比性高，讓處理后期的異常解釋更為直觀可靠。

這里利用Matlab,以聯(lián)剖法數(shù)據(jù)處理為例，通過(guò)編寫程序基本完成了聯(lián)剖法數(shù)據(jù)的全套處理，并且利用圖形用戶界面GUI編寫了一套可視化的聯(lián)剖法數(shù)據(jù)處理軟件，使處理過(guò)程更為簡(jiǎn)潔直觀[1]。也希望以此為窗口，啟發(fā)更多地質(zhì)工作者將Matlab運(yùn)用到地質(zhì)領(lǐng)域的各個(gè)方面中去。

1 聯(lián)剖法中的地形校正

1.1 地形校正計(jì)算方法

在電法勘探中,起伏的地形往往使有用異常面目皆非。因而在對(duì)電法數(shù)據(jù)資料進(jìn)行地質(zhì)解釋時(shí), 克服地形影響一直是重要問(wèn)題之一。為克服地形影響，需要計(jì)算求出由起伏地形引起的純地形影響 。這里利用Matlab,通過(guò)數(shù)值模擬來(lái)計(jì)算純地形影響。

選取的地改方法較為簡(jiǎn)單。即利用地表與水平面之間的夾角α，以及供電電極與記錄點(diǎn)的高程數(shù)據(jù)，根據(jù)視電阻率的微分表達(dá)式(式(1))[2]，計(jì)算求取純地形影響。

(1)

jMN*H′=jO*H

(2)

由此可得

(3)

進(jìn)而求得歸一化的純地形異常式(4)。

(4)

利用純地形異常即可根據(jù)式(4)求得進(jìn)過(guò)地改后的視電阻率值：

(5)

式中：ρS為野外實(shí)測(cè)電阻率；AO′、sin∠O′AO可根據(jù)供電電極、記錄點(diǎn)的位置及高程求得。

圖1 地形改正計(jì)算模型Fig.1 Calculation model for topographic correction

1.2 利用Matlab對(duì)不同地形進(jìn)行模擬計(jì)算的實(shí)例

圖2是用純地形影響求解方法求取的二維山脊和二維山谷地形上的聯(lián)合剖面視電阻率曲線。

圖2 模擬計(jì)算模型及純地形影響剖面圖Fig.2 Relief model and the result of pure topographic influence(a)山脊地形上聯(lián)剖剖面曲線； (b)山谷地形上聯(lián)剖剖面曲線

2 單線剖面圖繪制

在分析單條地電勘探線數(shù)據(jù)時(shí)，需要先繪制單線剖面圖。繪制剖面圖時(shí)可選取不同電性參數(shù)，從而從不同角度對(duì)地電測(cè)量結(jié)果進(jìn)行分析。這其中能反映測(cè)線電性特征的參數(shù)主要有以下幾種：①視電阻率；②視極化率；③比值參數(shù)；④岐離帶參數(shù)等。

在作者編寫的單線剖面圖Matlab程序[4]中，可以選擇是否進(jìn)行地形校正，并可以選取需要輸出的參數(shù)類型或者輸出全部參數(shù)曲線。下面根據(jù)實(shí)測(cè)數(shù)據(jù)給出具體實(shí)例。

2.1 地形校正

在做剖面圖時(shí)，當(dāng)?shù)匦纹鸱容^大時(shí)，會(huì)對(duì)測(cè)量結(jié)果產(chǎn)生較大影響，所以在繪制圖形前需要確定是否需要進(jìn)行地改處理[5]。

HIGH=xlsread('data.xlsx'); %讀取工區(qū)高程數(shù)據(jù)

circle=(AO-step/2)/step; %求AO間相隔的點(diǎn)數(shù)

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;

遠(yuǎn)程脈沖水表的動(dòng)力裝置是葉輪，當(dāng)水流通過(guò)葉輪時(shí)，帶動(dòng)葉輪旋轉(zhuǎn)，而水流的流速與葉輪的轉(zhuǎn)速成正比，因水流驅(qū)動(dòng)葉輪出噴口的截面積為常數(shù)，故葉輪的轉(zhuǎn)速與流量也成正比。因生產(chǎn)用水中含有少量雜質(zhì)，長(zhǎng)時(shí)間運(yùn)行后，水表葉輪會(huì)逐步聚集纏繞部分纖維物，導(dǎo)致水表不能正常工作。

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中，觀察地形可以發(fā)現(xiàn)，水平距離200 m的高差達(dá)到了50 m，這將對(duì)測(cè)量結(jié)果造成影響。

對(duì)比圖3(a),圖3(b)可以發(fā)現(xiàn)，曲線發(fā)生了變化。①?gòu)膱D3(a)中可以看到，視電阻率曲線245 m處存在一個(gè)低阻正交點(diǎn)，而從圖3(b)中可以看到，這一交點(diǎn)其實(shí)是不存在的，而是受地形起伏影響而產(chǎn)生的“假交點(diǎn)”。②圖3(a)中，410 m位置處存在一個(gè)低阻正交點(diǎn)，在圖3(b)中，這一交點(diǎn)出現(xiàn)在370 m處，這說(shuō)明地形起伏的影響使得交點(diǎn)的位置發(fā)生了變化。

圖3 地形校正曲線圖Fig.3 Profile of topographic correction(a)野外實(shí)測(cè)視電阻率剖面圖；(b)地改后的視電阻率剖面圖；(c)純地形影響剖面圖；(d)地形剖面圖

由圖3可以看出，當(dāng)?shù)匦纹鸱^大時(shí)，倘若不對(duì)原始數(shù)據(jù)進(jìn)行地形校正，直接根據(jù)交點(diǎn)性質(zhì)與位置加以解釋的話，會(huì)得到錯(cuò)誤的推斷。

2.2 比值參數(shù)

當(dāng)觀測(cè)結(jié)果受地表局部不均勻體影響時(shí)，觀測(cè)到的聯(lián)合剖面曲線將出現(xiàn)同步跳躍現(xiàn)象。這將會(huì)給基于交點(diǎn)性質(zhì)來(lái)解釋觀測(cè)結(jié)果的聯(lián)合剖面法造成一定的困難。為此，可以通過(guò)計(jì)算比值參數(shù)來(lái)對(duì)地表不均勻體的影響進(jìn)行壓制。比值參數(shù)按式(6)進(jìn)行計(jì)算。

(6)

求取比值參數(shù)和岐離帶參數(shù)的程序代碼如下：

Fa=psAC./psBC; %計(jì)算求取Fa

Fb=psBC./psAC; %計(jì)算求取Fb

subplot(2,1,1); %分幅繪制比值曲線

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

hold on;

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

for i=1:1(size(psAC)-1); %計(jì)算岐離帶參數(shù)

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

end

subplot(2,1,2); %分幅繪制岐離帶參數(shù)曲線

plot(x,lamuta,'r');

觀察圖4可以看出，與經(jīng)過(guò)地形改正后的視電阻率曲線進(jìn)行比較，雖然地形改正消除了曲線受地形的影響，但由于局部不均勻體的影響，視電阻率曲線靠的更近，同時(shí)同步跳躍。計(jì)算比值參數(shù)后，比值曲線分離較大，異常的交點(diǎn)更加明顯。

2.3 歧離帶參數(shù)

歧離帶參數(shù)是一種用來(lái)突出地質(zhì)異常的比值參數(shù)。對(duì)于低阻直立薄板上，聯(lián)合剖面的歧離帶參數(shù)曲線為單峰異常。模型模擬結(jié)果表明，該參數(shù)具有一定的復(fù)合異常分析能力。式(7)為岐離帶參數(shù)的計(jì)算公式。

(7)

圖5是利用野外實(shí)測(cè)數(shù)據(jù)通過(guò)計(jì)算得到的岐離帶參數(shù)曲線。在視電阻率曲線中，低阻正交點(diǎn)往往指示著低阻體的存在，而對(duì)比觀察岐離帶曲線與視電阻率曲線可以看出，剖面290 m、435 m、700 m、920 m出現(xiàn)的低阻正交點(diǎn)都有相應(yīng)的岐離帶單峰異常出現(xiàn)，這說(shuō)明岐離帶參數(shù)對(duì)于低阻體的存在確實(shí)有著指示作用。但需要注意的是，這一參數(shù)受干擾影響大，在實(shí)際運(yùn)用中需謹(jǐn)慎使用。

圖4 比值參數(shù)曲線圖Fig.4 Profile of ratio parameters(a)視電阻率曲線；(b)比值參數(shù)曲線

圖5 岐離帶參數(shù)曲線圖Fig.5 Profile of straggling parameters(a)視電阻率曲線；(b)岐離帶參數(shù)曲線

3 平面剖面圖繪制

平面剖面圖是直流電法勘探中一種可以同時(shí)顯示某種地球物理場(chǎng)剖面及其平面特征的重要圖件。

平剖圖中，各測(cè)線異常都用剖面形式表示，且各剖面比例尺統(tǒng)一。因而通過(guò)各測(cè)線之間的對(duì)比分析能清楚地表示異常平面和剖面的分布規(guī)律。進(jìn)一步探明礦區(qū)異常的走向及其延展性。

在作者編寫的Matlab平剖圖繪制程序中，可自行選擇需要輸出的剖面號(hào)，選擇是否進(jìn)行地形校正，并且可以選擇需要輸出的參數(shù)類型(包括視電阻率、視極化率、比值參數(shù)、岐離帶參數(shù))。程序輸出圖形，見(jiàn)圖6。對(duì)于這些輸出參數(shù)的控制，具體的Matab程序代碼如下：

%以矩陣形式輸出所需處理的測(cè)線號(hào)

%選擇是否地改并選擇需要輸出的曲線類型

inputlines=input('請(qǐng)輸入處理線號(hào)：');

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

CurveType=input('輸出曲線類型(視電阻率：1 比值曲線：2 岐離帶曲線：3)')

%根據(jù)輸入的線號(hào)編程循環(huán)輸出各剖面曲線

圖6 不同參數(shù)的平剖圖曲線Fig.6 Profile of different parameters(a)帶地改的視電阻率曲線；(b)不帶地改的比值參數(shù)曲線

4 GUI設(shè)計(jì)的可視化程序

為使數(shù)據(jù)處理及成圖過(guò)程更為直觀、便捷，利用Matlab GUI設(shè)計(jì)了以下可視化程序。該程序共有三個(gè)處理模塊：①單線剖面圖繪制；②平剖圖繪制；③點(diǎn)位投影。圖7為編寫的主菜單界面，用以選擇調(diào)用處理模塊。

圖7 聯(lián)剖法數(shù)據(jù)處理程序主菜單Fig.7 Main menu of the data processing software for composite profiling method

單線剖面繪制程序用以輸出單條地電斷面的剖面圖。程序首先需要選擇輸入Excel格式的數(shù)據(jù)文件及輸出的線號(hào)，即可輸出剖面圖。做圖過(guò)程中可以隨時(shí)更改輸出線號(hào)，輸出參數(shù)，是否地改，輸出剖面的范圍等。此外還可以設(shè)置圖表輸出參數(shù)來(lái)使圖像美觀。程序界面見(jiàn)圖8。

平剖圖繪制程序用以繪制平面剖面圖。選擇Excel格式的數(shù)據(jù)文件，并輸入需要輸出的所有測(cè)線號(hào)，即可繪制出平剖圖。同時(shí)可以選擇輸出的電性參數(shù)，是否地改，圖像會(huì)隨設(shè)定參數(shù)的改變隨時(shí)更新。其次還能設(shè)置圖表參數(shù)。程序界面見(jiàn)圖9。

點(diǎn)位投影用來(lái)輸出測(cè)區(qū)的三維地形，并利用視電阻率聯(lián)剖數(shù)據(jù)計(jì)算求得對(duì)稱四極裝置的視電阻率異常，將其投影到三維地形上并上色。其次還將投影所有測(cè)點(diǎn)位置，以及程序拾取得到的正反交點(diǎn)位置到三維地形上，以便更為直觀的分析觀察測(cè)量結(jié)果。程序界面見(jiàn)圖10。

圖8 單線剖面圖繪制程序Fig.8 Mapping program for single profile

圖9 平剖圖繪制程序Fig.9 Mapping program for profile map

5 結(jié)束語(yǔ)

通過(guò)Matlab編程基本實(shí)現(xiàn)了聯(lián)剖法數(shù)據(jù)整理及成圖的全套處理。而GUI編寫的可視化程序使得處理過(guò)程簡(jiǎn)單直觀。當(dāng)然所編寫的程序中地形改正方法還有待進(jìn)一步改善，數(shù)據(jù)文件的兼容性問(wèn)題也需要進(jìn)一步提高?？傊甅atlab強(qiáng)大的數(shù)據(jù)處理功能使得它在整個(gè)地學(xué)領(lǐng)域都有其用武之處。希望能看到Matlab在地質(zhì)工作中得到更為廣泛地使用。

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

[2] 傅良魁. 應(yīng)用地球物理教程—電法 放射性 地?zé)?[M]. 北京：地質(zhì)出版社，1991. FU L K. Applied geophysics [M]. Beijing : Geological Publishing House, 1991. (In Chinese)

圖10 點(diǎn)位投影程序Fig.10 Mapping program for point pojection

[3] 李金銘. 地電場(chǎng)與電法勘探 [M]. 北京：地質(zhì)出版社, 2005． LI J M. Electric and electrical prospecting [M].Beijing:Geological Publishing House,2005. (In Chinese)

[4] 童孝忠,王濤,柳建新. Matlab軟件包在地球物理數(shù)值計(jì)算與輔助教學(xué)中的應(yīng)用[J].物探化探計(jì)算技術(shù), 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] 湯井田,辛?xí)?王冉. 點(diǎn)電源下復(fù)雜角域地形影響及校正[J]. 吉林大學(xué)學(xué)報(bào)(地球科學(xué)版),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] 許瑞華. 比值參數(shù)在激電聯(lián)剖視極化率數(shù)據(jù)處理中的應(yīng)用[J]. 四川地質(zhì)學(xué)報(bào),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的工程物探軟件快速開(kāi)發(fā)[J]. 地球物理學(xué)進(jìn)展,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數(shù)據(jù)體三維可視化技術(shù)[J].地球物理學(xué)進(jìn)展,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

中國(guó)地質(zhì)調(diào)查項(xiàng)目(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