地质方向数据玫瑰花图自动绘制及其应用

余继峰，李政宏，陈 曦

(山东科技大学 地球科学与工程学院，山东 青岛 266590)

地质方向数据玫瑰花图自动绘制及其应用

余继峰，李政宏，陈 曦

(山东科技大学 地球科学与工程学院，山东 青岛 266590)

玫瑰花图是方向统计数据的主要表达方式之一，能直观地反映诸如节理、断层、沉积构造等定向地质数据的优选方向和总体分布特征，有助于分析解决区域构造应力场方向以及沉积盆地水流或物源方向等地质问题。过去大多采用人工统计手工绘制的方式实现，工作繁琐，效率低。计算机技术广泛应用以来，不少学者尝试使用不同的计算机语言编写程序解决绘图自动化的问题，但依然存在一些不足之处。本研究借助第四代编程语言MATLAB强大的数据处理功能和丰富的函数库，编写少量代码实现对该类地质数据的自动统计、玫瑰花图自动绘图和自动输出等功能。与用MATLAB自带的rose()函数绘制的玫瑰花图相比，更符合地质工作者的使用习惯，通过鄂东临汾地区燕山期节理玫瑰花图的自动绘制，展示了该方法的实用性。

玫瑰花图；方向数据；MATLAB语言；优选方位

玫瑰花图在地质领域中应用广泛，清晰形象。但人工绘制玫瑰花图不仅费时费力，统计量大，而且容易出错，成图不易保存。因此，在计算机尚未普及时便有学者用PC-1500袖珍计算机借助BASIC算法语言[1]进行玫瑰花图自动绘制方面的研究，后来又有学者编写了玫瑰花图自动绘制软件GSMROSE[2]。 Windows操作系统出现后，陆续出现了用 Visual Basic 语言开发的古水流玫瑰花图绘制软件PC99[3]。近年来，又陆续出现了基于CAD[4]、EXCEL[5]和面向对象的程序设计语言VB[6-7]、VC++[8]编写的玫瑰花图绘制程序。然而，上述方法一方面存在编程语言复杂、程序需要从底层开发、出错率高、数据储存量少等问题，另一方面图形存储格式受限、使用不便。

被誉为第四代编程语言的MATLAB出现后，因其工程计算库函数丰富等特点为工程领域数值计算与数据处理以及地质领域研究提供了一种新方法[9-12]。另外，其程序语言是一种基于矩阵的高级语言，具有以下几个特点：功能强大的数值运算功能、强大的图形处理能力、高级但简单的程序环境、丰富的工具箱等，只需少量代码就能实现其他语言较难实现的多种功能。张立波[13]曾应用MATLAB编程绘制风向频率玫瑰图，而在地质领域方向数据分析[14-15]方面的应用并不多见。应用MATLAB语言编写了玫瑰花图绘制程序，并通过实例数据验证了其良好的实用性，提高了地学科技人员的工作效率。

1 程序处理流程和设计步骤

运用MATLAB设计玫瑰花图程序，主要由dip()、circle()、strike()、semicircle()四个函数组成。circle()和semicircle()分别用于绘制倾向和走向玫瑰花图的轮廓。dip()用于统计倾向数据并绘制倾向玫瑰花图，strike()用于统计走向数据并绘制走向玫瑰花图。

步骤一：加载并处理数据

1) 加载数据：

源数据文件要求使用文本文件,数据格式如下：

10,145,270,353……

从已保存的文件directional_1.txt中加载数据的语句如下：

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

forj=1:length(varargin)

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

end%此循环用于统计输入数据的数目，并将原始数据保存在data变量中

direction=0:pi/18:35*pi/18;%初始化36个方向及其对应弧度

r1=0.7; font_size=7;% r1为标记文本与外轮廓的距离，font_size为玫瑰花图中标记文本字号

2) 处理数据

i.倾向范围在0°～360°，玫瑰花图用整圆，用原始数据即可。

ii.走向玫瑰花图用上半圆，需要将大于90°小于等于270°的数据调整在上半圆范围内，处理代码如下：

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

data(find((data>180)&(data<=270)))=[];%将被修改的数据清除

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

data(find((data>90)&(data<=180)))=[];%将被修改的数据清除

newdata=cat(2,data,newdata1,newdata2);%修改后的数据保存在新变量中

步骤二：将原始数据按花瓣大小分组统计，统计结果以表格形式保存

1) 倾向数据统计代码

head1=1;head2=10;%head1为第一个分组倾向的最小值，head2为最大值

fori=1:36 %以10°为一个分组，共分36组

n(i)=sum(((data>=head1)&(data<=head2)),2);%统计每个分组的数据数目

ave(i)=round(mean(data(find((data>=head1)&(data<=head2)))));%计算每个分组内数据的平均值

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

end%此循环用于计算每个分组内数据数目和平均方向

调用xlswrite函数把统计好的数据写入Excel文件，调用格式：

xlswrite(filename,M,range);%filename为目标文件名，M为写入Excel的数据矩阵(n(i)、ave(i))，range 为写入的单元格区域(n(i)写入B2至B37,ave(i)写入C2至C37)

2) 走向数据统计代码

走向只需要分别统计1°～90°和271°～360°内的数据即可，不再赘述。

步骤三：绘制玫瑰花图

1) 需要分别绘制轮廓线、方位线、方位角度，依据自动统计的数据绘制玫瑰花瓣。以绘制倾向玫瑰花图为例：

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

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

num=cat(2,newnum1′,newnum2′);%提取各分组内对应数据数目和平均方向，并将列数据变形成行数据

nmax=max(n);%统计出n(i)中数据最多的那一组数据的数目

2) 设置玫瑰花图属性(字体大小、颜色、位置)

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) 绘制从原点发出的36条方向线

forii=1:numel(direction)

hold on;

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

end

5) 在圆上标记方位

由于各个方位角度在圆上标注的位置不同，因此标注文本的水平对齐方式应不同。以标记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,′节理倾向玫瑰花图′,′jpg′);

2 应用实例

2.1 自编程序作图

应用上述程序，根据文献[16]中鄂东临汾地区燕山期共轭剪节理实测数据(表1)，在MATLAB软件平台上分别绘制了节理走向、倾向玫瑰花图(图1、图2)。该图清晰地展示出燕山中期共轭剪节理系走向主要为NEE和NW向，及近NS和EW向展布，平均优选方位分别为55°和325°，355°和85°，指示研究区遭受过NWW向及NW向挤压应力作用。由于研究区张节理较少，地层平缓(<15°)构造简单，普遍发育高角度剪节理，为早期平面共轭剪节理[17]，结合共轭剪节理相互交切形态，恢复燕山中期构造应力以NW向为主，局部转为NWW向，燕山晚期为NWW向。

图1 燕山期共轭剪节理走向玫瑰花图Fig.1 Strike rose diagram of Yanshanian conjugated shear joints

图2 燕山期共轭剪节理倾向玫瑰花图Fig.2 Dip rose diagram of Yanshanian conjugated shear joints表1 燕山期共轭剪节理原始数据Tab.1 Original data Yanshan epoch conjugated shear joints

观测点区位点号共轭剪节理倾向/(°)倾角/(°)倾向/(°)倾角/(°)形成期次(挤压应力方向)观测点区位点号共轭剪节理倾向/(°)倾角/(°)倾向/(°)倾角/(°)形成期次(挤压应力方向)东南部17428627787晚期(NNW)西部73118410181中期(NW)东南部141518623486中期(NWW)西部49908818187中期(NW)东南部162368332882中期(NWW)西部70189769079中期(NW)东南部4307231085中期(NWW)西部531048410783中期(NW)东南部11708623985中期(NW)窑渠背斜1292257913284中期(NWW)

续表

数据统计结果自动保存在Excel文件中(表2)，方便查看。

表2 燕山期共轭剪节理统计数据

续表

2.2 MATLAB库函数中rose()作图

MATLAB自带的函数rose()也可以进行玫瑰花图的制作[18]。函数rose使用逆时针方向计数，沿着直角坐标x轴正方向从0°开始。但在地球科学中，0°点表示北，90°点表示东，并且角度按照顺时针方向增加。命令view通过+90°(方位角)旋转图形，通过-90°(仰角)镜像图形[19-20]。用上述实例数据重新绘制倾向、走向玫瑰花图(图3、图4)。

图3 燕山期共轭剪节理走向玫瑰花图Fig.3 Strike rose diagram of Yanshanian conjugated shear joints

图4 燕山期共轭剪节理倾向玫瑰花图Fig.4 Dip rose diagram of Yanshanian conjugated shear joints

图5 燕山期共轭剪节理倾向玫瑰花图Fig.5 Dip rose diagram of Yanshanian conjugated shear joints

具体代码如下：

clear

data_degrees=load(‘directional_1.txt’);

data_radians=pi*data_degress/180;%将角度值转换成弧度值

rose(data_radians,24);

view(90.-90)

通过对比可以发现，rose()函数绘制的玫瑰花图与原始玫瑰花图相比优势方位不明显，虽然编写的代码数量少，对走向数据进行成图时需要人工校正在1°～90°和271°～360°范围内，增加了工作量和出错率。再者rose()函数不具有将数据自动统计输出的特点，而本程序实现了将数据自动统计并自动保存成Excel文件。

2.3 不同花瓣宽度的比较

对于实例数据倾向玫瑰花图来说，选用15°间隔作图(图1)效果好于10°间隔作图(图5)，前者能清晰展示出研究区发育的两组共轭剪节理方向，即NEE向和NW向，及近NS和EW向，而后者由于花瓣数量过多，在一定程度上影响了对区域应力场的正确判断分析。

3 结语

1) 借助MATLAB强大的数据处理功能和丰富的函数库，编写少量代码实现了对该类地质方向数据的自动统计、玫瑰花图自动绘图和自动输出等功能。但仅用rose()库函数制作玫瑰花图除有时影响优选方位判断外，不符合地质工作者对玫瑰花图的使用习惯。

2) 如需对玫瑰花图填色，只要利用fill函数对不同方位充填自己满意的颜色即可。如需输出其他格式图形，只要设置saveas函数的’format’属性即可。本文保存的是jpg格式。

3) 玫瑰花图花瓣的半径直接表示统计结果的多寡，体现了优选方向。花瓣宽度间隔可在5°～30°之间选择，玫瑰花图的花瓣宽度越大，显示结果的分辨率和复杂程度就越低，宽度小时更能表现出细节，但复杂程度增高，有时会影响对优选方位的判断，一般间隔选用10°～15°，可视具体情况而定。

[1]文朴.用PC-1500袖珍电子计算机统计节理并绘制玫瑰花图[J].地质与勘探,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].华东地质学院学报,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.图形显示和比较古水流数据的一种软件( PC99):以青藏高原北部可可西里盆地新生代古水流数据为例[J].沉积学报,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]段其所.节理裂隙玫瑰花图程序设计及在水电工程中的应用[J].云南水力发电,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 绘制节理走向玫瑰花图[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]张俊林,杜丙义,李建民,等.应用VB编程实现气象资料的自动统计及风向玫瑰图的绘制[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]况源,周小明,梁富强.基于VC++的风玫瑰图绘制程序设计与实现[J].电子设计工程,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]余继峰,于泳,付文钊,等.测井数据Matlab插值与地质旋回性分析应用[J].煤炭学报,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的测井曲线频谱分析及其在地质研究中的应用：以川东北地区二叠系长兴组为例[J].天然气地球科学,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的小波变换在沉积旋回研究中的应用[J].物探化探计算技术,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的风玫瑰图绘制[J].电脑编程技巧与维护,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]王琳琳.煤储层节理发育非均质性评价的构造动力学方法及其应用：以鄂东临汾地区为例[D].徐州:中国矿业大学,2014:53-59.

[17]朱昊.鄂尔多斯盆地西缘中南段地质结构及其形成演化[D].北京：中国地质大学,2015:63-66.

[18] TRAUTH .M著，宋久旭，等 译.地球科学中的MATLAB应用[M].北京：国防工业出版社，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.

(责任编辑：高丽华)

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

国家自然科学基金项目(41472092)；国家级大学生创新创业训练计划项目(201510424005,201510424003)

余继峰(1964—)，男，安徽萧县人，教授，博士生导师，主要从事能源盆地分析和测井地质定量解释、旋回地层学等方面的研究.E-mail:yujifeng05@163.com 李政宏(1994—)，男，山东青岛人，硕士研究生，主要从事计算机技术在地质领域应用方面的研究,本文通信作者. E-mail:2446347999@qq.com

P623.6

A

1672-3767(2017)03-0001-08