地質(zhì)方向數(shù)據(jù)玫瑰花圖自動繪制及其應(yīng)用

余繼峰，李政宏，陳 曦

(山東科技大學(xué) 地球科學(xué)與工程學(xué)院，山東 青島 266590)

地質(zhì)方向數(shù)據(jù)玫瑰花圖自動繪制及其應(yīng)用

余繼峰，李政宏，陳 曦

(山東科技大學(xué) 地球科學(xué)與工程學(xué)院，山東 青島 266590)

玫瑰花圖是方向統(tǒng)計數(shù)據(jù)的主要表達(dá)方式之一，能直觀地反映諸如節(jié)理、斷層、沉積構(gòu)造等定向地質(zhì)數(shù)據(jù)的優(yōu)選方向和總體分布特征，有助于分析解決區(qū)域構(gòu)造應(yīng)力場方向以及沉積盆地水流或物源方向等地質(zhì)問題。過去大多采用人工統(tǒng)計手工繪制的方式實現(xiàn)，工作繁瑣，效率低。計算機(jī)技術(shù)廣泛應(yīng)用以來，不少學(xué)者嘗試使用不同的計算機(jī)語言編寫程序解決繪圖自動化的問題，但依然存在一些不足之處。本研究借助第四代編程語言MATLAB強(qiáng)大的數(shù)據(jù)處理功能和豐富的函數(shù)庫，編寫少量代碼實現(xiàn)對該類地質(zhì)數(shù)據(jù)的自動統(tǒng)計、玫瑰花圖自動繪圖和自動輸出等功能。與用MATLAB自帶的rose()函數(shù)繪制的玫瑰花圖相比，更符合地質(zhì)工作者的使用習(xí)慣，通過鄂東臨汾地區(qū)燕山期節(jié)理玫瑰花圖的自動繪制，展示了該方法的實用性。

玫瑰花圖；方向數(shù)據(jù)；MATLAB語言；優(yōu)選方位

玫瑰花圖在地質(zhì)領(lǐng)域中應(yīng)用廣泛，清晰形象。但人工繪制玫瑰花圖不僅費時費力，統(tǒng)計量大，而且容易出錯，成圖不易保存。因此，在計算機(jī)尚未普及時便有學(xué)者用PC-1500袖珍計算機(jī)借助BASIC算法語言[1]進(jìn)行玫瑰花圖自動繪制方面的研究，后來又有學(xué)者編寫了玫瑰花圖自動繪制軟件GSMROSE[2]。 Windows操作系統(tǒng)出現(xiàn)后，陸續(xù)出現(xiàn)了用 Visual Basic 語言開發(fā)的古水流玫瑰花圖繪制軟件PC99[3]。近年來，又陸續(xù)出現(xiàn)了基于CAD[4]、EXCEL[5]和面向?qū)ο蟮某绦蛟O(shè)計語言VB[6-7]、VC++[8]編寫的玫瑰花圖繪制程序。然而，上述方法一方面存在編程語言復(fù)雜、程序需要從底層開發(fā)、出錯率高、數(shù)據(jù)儲存量少等問題，另一方面圖形存儲格式受限、使用不便。

被譽為第四代編程語言的MATLAB出現(xiàn)后，因其工程計算庫函數(shù)豐富等特點為工程領(lǐng)域數(shù)值計算與數(shù)據(jù)處理以及地質(zhì)領(lǐng)域研究提供了一種新方法[9-12]。另外，其程序語言是一種基于矩陣的高級語言，具有以下幾個特點：功能強(qiáng)大的數(shù)值運算功能、強(qiáng)大的圖形處理能力、高級但簡單的程序環(huán)境、豐富的工具箱等，只需少量代碼就能實現(xiàn)其他語言較難實現(xiàn)的多種功能。張立波[13]曾應(yīng)用MATLAB編程繪制風(fēng)向頻率玫瑰圖，而在地質(zhì)領(lǐng)域方向數(shù)據(jù)分析[14-15]方面的應(yīng)用并不多見。應(yīng)用MATLAB語言編寫了玫瑰花圖繪制程序，并通過實例數(shù)據(jù)驗證了其良好的實用性，提高了地學(xué)科技人員的工作效率。

1 程序處理流程和設(shè)計步驟

運用MATLAB設(shè)計玫瑰花圖程序，主要由dip()、circle()、strike()、semicircle()四個函數(shù)組成。circle()和semicircle()分別用于繪制傾向和走向玫瑰花圖的輪廓。dip()用于統(tǒng)計傾向數(shù)據(jù)并繪制傾向玫瑰花圖，strike()用于統(tǒng)計走向數(shù)據(jù)并繪制走向玫瑰花圖。

步驟一：加載并處理數(shù)據(jù)

1) 加載數(shù)據(jù)：

源數(shù)據(jù)文件要求使用文本文件,數(shù)據(jù)格式如下：

10,145,270,353……

從已保存的文件directional_1.txt中加載數(shù)據(jù)的語句如下：

clear;data=load(‘directional_1.txt’);

forj=1:length(varargin)

data(j)=varargin{j}(:);

end%此循環(huán)用于統(tǒng)計輸入數(shù)據(jù)的數(shù)目，并將原始數(shù)據(jù)保存在data變量中

direction=0:pi/18:35*pi/18;%初始化36個方向及其對應(yīng)弧度

r1=0.7; font_size=7;% r1為標(biāo)記文本與外輪廓的距離，font_size為玫瑰花圖中標(biāo)記文本字號

2) 處理數(shù)據(jù)

i.傾向范圍在0°～360°，玫瑰花圖用整圓，用原始數(shù)據(jù)即可。

ii.走向玫瑰花圖用上半圓，需要將大于90°小于等于270°的數(shù)據(jù)調(diào)整在上半圓范圍內(nèi)，處理代碼如下：

data1=data(find((data>180)&(data<=270)));newdata1=data1-180;

data(find((data>180)&(data<=270)))=[];%將被修改的數(shù)據(jù)清除

data2=data(find((data>90)&(data<=180)));newdata2=data2+180;

data(find((data>90)&(data<=180)))=[];%將被修改的數(shù)據(jù)清除

newdata=cat(2,data,newdata1,newdata2);%修改后的數(shù)據(jù)保存在新變量中

步驟二：將原始數(shù)據(jù)按花瓣大小分組統(tǒng)計，統(tǒng)計結(jié)果以表格形式保存

1) 傾向數(shù)據(jù)統(tǒng)計代碼

head1=1;head2=10;%head1為第一個分組傾向的最小值，head2為最大值

fori=1:36 %以10°為一個分組，共分36組

n(i)=sum(((data>=head1)&(data<=head2)),2);%統(tǒng)計每個分組的數(shù)據(jù)數(shù)目

ave(i)=round(mean(data(find((data>=head1)&(data<=head2)))));%計算每個分組內(nèi)數(shù)據(jù)的平均值

head1=head1+10;head2=head2+10;

end%此循環(huán)用于計算每個分組內(nèi)數(shù)據(jù)數(shù)目和平均方向

調(diào)用xlswrite函數(shù)把統(tǒng)計好的數(shù)據(jù)寫入Excel文件，調(diào)用格式：

xlswrite(filename,M,range);%filename為目標(biāo)文件名，M為寫入Excel的數(shù)據(jù)矩陣(n(i)、ave(i))，range 為寫入的單元格區(qū)域(n(i)寫入B2至B37,ave(i)寫入C2至C37)

2) 走向數(shù)據(jù)統(tǒng)計代碼

走向只需要分別統(tǒng)計1°～90°和271°～360°內(nèi)的數(shù)據(jù)即可，不再贅述。

步驟三：繪制玫瑰花圖

1) 需要分別繪制輪廓線、方位線、方位角度，依據(jù)自動統(tǒng)計的數(shù)據(jù)繪制玫瑰花瓣。以繪制傾向玫瑰花圖為例：

num1=B(2:10);newnum1=flipud(num1);

num2=B(11:37);newnum2=flipud(num2);

num=cat(2,newnum1′,newnum2′);%提取各分組內(nèi)對應(yīng)數(shù)據(jù)數(shù)目和平均方向，并將列數(shù)據(jù)變形成行數(shù)據(jù)

nmax=max(n);%統(tǒng)計出n(i)中數(shù)據(jù)最多的那一組數(shù)據(jù)的數(shù)目

2) 設(shè)置玫瑰花圖屬性(字體大小、顏色、位置)

f1=figure(2);

set(f1,′units′,′normalized′,′position′,[0.10,0.10,0.80,0.80],′color′,′w′);

axis equal; box off; set(gca,′visible′,′off′); hold on;

3) 畫圓

r=nmax*5; circle(r);

4) 繪制從原點發(fā)出的36條方向線

forii=1:numel(direction)

hold on;

plot([0,r*cos(direction(ii))],[0,r*sin(direction(ii))],′k:′,′LineWidth′,0.5);

end

5) 在圓上標(biāo)記方位

由于各個方位角度在圓上標(biāo)注的位置不同，因此標(biāo)注文本的水平對齊方式應(yīng)不同。以標(biāo)記0°～90°為例：

forii=1:9

text((r+r1)*cos(direction(ii)),(r+r1)*sin(direction(ii)),strcat(fangxiang(ii)),’fontsize’,font_size,′horizontalAlignment′,′left′);

6) 繪制玫瑰花瓣

forii=1:numel(direction)

if(ii

hold on;

plot([num(ii)*5*cos(direction(ii)+pi/36),num(ii+1)*5*cos(direction(ii+1)+pi/36)],[num(ii)*5*sin(direction(ii)+pi/36),num(ii+1)*5*sin(direction(ii+1)+pi/36)],′k-′,′LineWidth′,1.5,′color′,′b′);

else

plot([num(36)*5*cos(direction(36)+pi/36),num(1)*5*cos(direction(1)+pi/36)],[num(36)*5*sin(direction(36)+pi/36),num(1)*5*sin(direction(1)+pi/36)],′k-′,′LineWidth′,1.5,′color′,′b′);

end

end

走向玫瑰花圖繪制過程與傾向玫瑰花圖基本一致，只需分別繪制1°～90°和271°～360°的花瓣即可，在此不再贅述。

步驟四：輸出圖形

set(f1,′inverthardcopy′,′off′); saveas(f1,′節(jié)理傾向玫瑰花圖′,′jpg′);

2 應(yīng)用實例

2.1 自編程序作圖

應(yīng)用上述程序，根據(jù)文獻(xiàn)[16]中鄂東臨汾地區(qū)燕山期共軛剪節(jié)理實測數(shù)據(jù)(表1)，在MATLAB軟件平臺上分別繪制了節(jié)理走向、傾向玫瑰花圖(圖1、圖2)。該圖清晰地展示出燕山中期共軛剪節(jié)理系走向主要為NEE和NW向，及近NS和EW向展布，平均優(yōu)選方位分別為55°和325°，355°和85°，指示研究區(qū)遭受過NWW向及NW向擠壓應(yīng)力作用。由于研究區(qū)張節(jié)理較少，地層平緩(<15°)構(gòu)造簡單，普遍發(fā)育高角度剪節(jié)理，為早期平面共軛剪節(jié)理[17]，結(jié)合共軛剪節(jié)理相互交切形態(tài)，恢復(fù)燕山中期構(gòu)造應(yīng)力以NW向為主，局部轉(zhuǎn)為NWW向，燕山晚期為NWW向。

圖1 燕山期共軛剪節(jié)理走向玫瑰花圖Fig.1 Strike rose diagram of Yanshanian conjugated shear joints

圖2 燕山期共軛剪節(jié)理傾向玫瑰花圖Fig.2 Dip rose diagram of Yanshanian conjugated shear joints表1 燕山期共軛剪節(jié)理原始數(shù)據(jù)Tab.1 Original data Yanshan epoch conjugated shear joints

觀測點區(qū)位點號共軛剪節(jié)理傾向/(°)傾角/(°)傾向/(°)傾角/(°)形成期次(擠壓應(yīng)力方向)觀測點區(qū)位點號共軛剪節(jié)理傾向/(°)傾角/(°)傾向/(°)傾角/(°)形成期次(擠壓應(yīng)力方向)東南部17428627787晚期(NNW)西部73118410181中期(NW)東南部141518623486中期(NWW)西部49908818187中期(NW)東南部162368332882中期(NWW)西部70189769079中期(NW)東南部4307231085中期(NWW)西部531048410783中期(NW)東南部11708623985中期(NW)窯渠背斜1292257913284中期(NWW)

續(xù)表

數(shù)據(jù)統(tǒng)計結(jié)果自動保存在Excel文件中(表2)，方便查看。

表2 燕山期共軛剪節(jié)理統(tǒng)計數(shù)據(jù)

續(xù)表

2.2 MATLAB庫函數(shù)中rose()作圖

MATLAB自帶的函數(shù)rose()也可以進(jìn)行玫瑰花圖的制作[18]。函數(shù)rose使用逆時針方向計數(shù)，沿著直角坐標(biāo)x軸正方向從0°開始。但在地球科學(xué)中，0°點表示北，90°點表示東，并且角度按照順時針方向增加。命令view通過+90°(方位角)旋轉(zhuǎn)圖形，通過-90°(仰角)鏡像圖形[19-20]。用上述實例數(shù)據(jù)重新繪制傾向、走向玫瑰花圖(圖3、圖4)。

圖3 燕山期共軛剪節(jié)理走向玫瑰花圖Fig.3 Strike rose diagram of Yanshanian conjugated shear joints

圖4 燕山期共軛剪節(jié)理傾向玫瑰花圖Fig.4 Dip rose diagram of Yanshanian conjugated shear joints

圖5 燕山期共軛剪節(jié)理傾向玫瑰花圖Fig.5 Dip rose diagram of Yanshanian conjugated shear joints

具體代碼如下：

clear

data_degrees=load(‘directional_1.txt’);

data_radians=pi*data_degress/180;%將角度值轉(zhuǎn)換成弧度值

rose(data_radians,24);

view(90.-90)

通過對比可以發(fā)現(xiàn)，rose()函數(shù)繪制的玫瑰花圖與原始玫瑰花圖相比優(yōu)勢方位不明顯，雖然編寫的代碼數(shù)量少，對走向數(shù)據(jù)進(jìn)行成圖時需要人工校正在1°～90°和271°～360°范圍內(nèi)，增加了工作量和出錯率。再者rose()函數(shù)不具有將數(shù)據(jù)自動統(tǒng)計輸出的特點，而本程序?qū)崿F(xiàn)了將數(shù)據(jù)自動統(tǒng)計并自動保存成Excel文件。

2.3 不同花瓣寬度的比較

對于實例數(shù)據(jù)傾向玫瑰花圖來說，選用15°間隔作圖(圖1)效果好于10°間隔作圖(圖5)，前者能清晰展示出研究區(qū)發(fā)育的兩組共軛剪節(jié)理方向，即NEE向和NW向，及近NS和EW向，而后者由于花瓣數(shù)量過多，在一定程度上影響了對區(qū)域應(yīng)力場的正確判斷分析。

3 結(jié)語

1) 借助MATLAB強(qiáng)大的數(shù)據(jù)處理功能和豐富的函數(shù)庫，編寫少量代碼實現(xiàn)了對該類地質(zhì)方向數(shù)據(jù)的自動統(tǒng)計、玫瑰花圖自動繪圖和自動輸出等功能。但僅用rose()庫函數(shù)制作玫瑰花圖除有時影響優(yōu)選方位判斷外，不符合地質(zhì)工作者對玫瑰花圖的使用習(xí)慣。

2) 如需對玫瑰花圖填色，只要利用fill函數(shù)對不同方位充填自己滿意的顏色即可。如需輸出其他格式圖形，只要設(shè)置saveas函數(shù)的’format’屬性即可。本文保存的是jpg格式。

3) 玫瑰花圖花瓣的半徑直接表示統(tǒng)計結(jié)果的多寡，體現(xiàn)了優(yōu)選方向。花瓣寬度間隔可在5°～30°之間選擇，玫瑰花圖的花瓣寬度越大，顯示結(jié)果的分辨率和復(fù)雜程度就越低，寬度小時更能表現(xiàn)出細(xì)節(jié)，但復(fù)雜程度增高，有時會影響對優(yōu)選方位的判斷，一般間隔選用10°～15°，可視具體情況而定。

[1]文樸.用PC-1500袖珍電子計算機(jī)統(tǒng)計節(jié)理并繪制玫瑰花圖[J].地質(zhì)與勘探,1984,28(11):40-43. WHEN Pu.Using PC-1500 pocket computer counts the joints and draws the rose diagram[J].Geology and Prospecting,1984,28(11):40-43.

[2]SELNER G,TAYLOR R.繪制玫瑰圖的程序[J].華東地質(zhì)學(xué)院學(xué)報,1991,14(3):285-288. SELNER G,TAYLOR R.Drawing the program of rose diagram[J].Journal of East China College of Geology,1991,14(3):285-288.

[3]劉志飛,LACHLAN K ,STEWART.圖形顯示和比較古水流數(shù)據(jù)的一種軟件( PC99):以青藏高原北部可可西里盆地新生代古水流數(shù)據(jù)為例[J].沉積學(xué)報,2002,20(2):354-358. LIU Zhifei,LACHLAN K,STEWART A.Software tool for graphically displaying and comparing paleocurrent data (PC99):An example utilizing paleocurrent data of the Cenozoic Hoh Xil basin,northern Tibet[J].Acta Sedimentologica Sinica,2002,20(2):354-358.

[4]段其所.節(jié)理裂隙玫瑰花圖程序設(shè)計及在水電工程中的應(yīng)用[J].云南水力發(fā)電,2012，28(4):56-57,146. DUAN Qisuo.Joint fissure rose diagram programming and the application in water-power engineering[J].Yunnan Water Power,2012,28(4):56-57,146.

[5]陳軍.基于 Excel 繪制節(jié)理走向玫瑰花圖[J].江淮水利科技,2014,36(6):10-11. CHEN Jun.Drawing joint strike rose diagram on the basis of Excel[J].Jianghuai Water Resource Science and Technology,2014,36(6):10-11.

[6]于國芳,閆寶珍,賁旭東,等.基于 VB 的古流向玫瑰花圖的繪制[J].西部探礦工程,2015,27(9):71-72. YU Guofang,YAN Baozhen,BEN Xudong,et al.Drawing paleocurrent rose digram on the basis of VB[J].West-China Exploration Engineering,2015,27(9):71-72.

[7]張俊林,杜丙義,李建民,等.應(yīng)用VB編程實現(xiàn)氣象資料的自動統(tǒng)計及風(fēng)向玫瑰圖的繪制[J].氣象水文海洋儀器,2014,31(3):63-65. ZHANG Junlin,DU Bingyi,LI Jianmin,et al.Application of VB in automatic statistics of meteorological data and the drawing of wind rose diagram[J].Meteorological,Hydrological and Marine Instruments,2014,31(3):63-65.

[8]況源,周小明,梁富強(qiáng).基于VC++的風(fēng)玫瑰圖繪制程序設(shè)計與實現(xiàn)[J].電子設(shè)計工程,2015,23(17):189-191. KUANG Yuan,ZHOU Xiaoming,LIANG Fuqiang.Based on VC++ wind rose drawing program design and implementation[J].Electronic Design Engineering,2015,23(17):189-191.

[9]余繼峰,于泳,付文釗,等.測井?dāng)?shù)據(jù)Matlab插值與地質(zhì)旋回性分析應(yīng)用[J].煤炭學(xué)報,2011,36(10):1679-1682. YU Jifeng,YU Yong,FU Wenzhao,et al.Application of interpolation of well logs based on Matlab to analysis of geological cyclicity[J].Journal of China Coal Society,2011,36(10):1679-1682.

[10]李鳳杰,趙俊興.基于Matlab的測井曲線頻譜分析及其在地質(zhì)研究中的應(yīng)用：以川東北地區(qū)二疊系長興組為例[J].天然氣地球科學(xué),2007,18(4):531-534. LI Fengjie,ZHAO Junxing.Spectrum analysis of logging curves using Matlab and its application in geology[J].Natural Gas Geoscience,2007,18(4):531-534.

[11]高迪,郭變青,邵龍義,等.基于Matlab的小波變換在沉積旋回研究中的應(yīng)用[J].物探化探計算技術(shù),2012,34(4):444-448. GAO Di,GUO Bianqing,SHAO Longyi,et al.Application of wavelet transform in sedimentary cycle research based on Matlab[J].Computing Techniques for Geophysical and Geochemical Exploration,2012,34(4):444-448.

[12]HOLZBERCHER E.Environmental modeling using Matlab[M].Berlin Heidelberg:Springer,2007:286-291.

[13]張立波.基于MATLAB的風(fēng)玫瑰圖繪制[J].電腦編程技巧與維護(hù),2012,19(18):26-27. ZHANG Libo.Plotting wind rose based on the Matlab program[J].Programming Skills and Maintenance,2012,19(18):26-27.

[14]BORRADAILE G.Statistics of earth science data:Their distribution in time,space and orientation[M].Berlin Heidelberg:Springer,2007:1-10.

[15]PRESS W H,TEUKOLSKY S A,VETTERLING W T,et al.Numerical recipes:The art of scientific computing 3rd edition[M].Cambrige:Cambrige University Press,2007:33-54.

[16]王琳琳.煤儲層節(jié)理發(fā)育非均質(zhì)性評價的構(gòu)造動力學(xué)方法及其應(yīng)用：以鄂東臨汾地區(qū)為例[D].徐州:中國礦業(yè)大學(xué),2014:53-59.

[17]朱昊.鄂爾多斯盆地西緣中南段地質(zhì)結(jié)構(gòu)及其形成演化[D].北京：中國地質(zhì)大學(xué),2015:63-66.

[18] TRAUTH .M著，宋久旭，等 譯.地球科學(xué)中的MATLAB應(yīng)用[M].北京：國防工業(yè)出版社，2015,7：120-180.

[19]PALM W J.Introduction to MATLAB7 for engineers[M].New York:McGraw-Hill,2004:264-267.

[20]DAVIS J C.Statistics and data analysis in geology[M].New York:John Wiley and Sons,2002:160-187.

(責(zé)任編輯：高麗華)

Automatic Plotting of Rose Diagram with Geological Data and Its Application

YU Jifeng,LI Zhenghong,CHEN Xi

( College of Earth Science and Engineering,Shandong University of Science and Technology, Qingdao,Shandong 266590,China )

Rose diagram,one of the major ways to visualize directional geological data,can clearly and intuitively reflect the dominant direction and general distribution of such geological data as the strike of joints,faults,and sedimentary structures.It is helpful in solving such geological problems as the direction of the regional tectonic stress field,the current or provenance directions of sedimentary basins.The rose diagram was commonly drawn manually,which is time consuming and with low efficiency.With the wide use of computer technology,different computer languages are used to design programs to solve the problem of automatic plotting,but there are still some disadvantages in many aspects.The authors of this paper wrote several lines of computer routines with the help of the ample function library and powerful data processing function of the fourth-generation programming language MATLAB platform to realize the automatic statistics of geological data,the automatic plotting of the rose diagram and the automatic output.Compared with the rose diagram plotted by the embedding function rose() of MATLAB,the automatically plotted rose program is more preferable to geologists.The automatic plotting of the rose diagram of Yanshanian joints demonstrates its validity and practicability.

rose diagram; directional data; MATLAB program; dominant direction

2016-09-14

國家自然科學(xué)基金項目(41472092)；國家級大學(xué)生創(chuàng)新創(chuàng)業(yè)訓(xùn)練計劃項目(201510424005,201510424003)

余繼峰(1964—)，男，安徽蕭縣人，教授，博士生導(dǎo)師，主要從事能源盆地分析和測井地質(zhì)定量解釋、旋回地層學(xué)等方面的研究.E-mail:yujifeng05@163.com 李政宏(1994—)，男，山東青島人，碩士研究生，主要從事計算機(jī)技術(shù)在地質(zhì)領(lǐng)域應(yīng)用方面的研究,本文通信作者. E-mail:2446347999@qq.com

P623.6

A

1672-3767(2017)03-0001-08