地质力学学报  2009, Vol. 15 Issue (3): 218-225
引用本文
李勃, 张廷山, 王占磊, 董银磊. 基于GIS的颗粒形状和应变分析[J]. 地质力学学报, 2009, 15(3): 218-225.
LI Bo, ZHANG Ting-shan, WANG Zhan-lei, DONG Yin-lei. ANALYSES ON GRAIN SHAPE AND LIMITED RESPONSE USING GIS[J]. Journal of Geomechanics, 2009, 15(3): 218-225.
基于GIS的颗粒形状和应变分析
李勃 , 张廷山 , 王占磊 , 董银磊     
西南石油大学资环学院, 四川 成都 610500
摘要:将GIS应用到岩石颗粒形状和有限应变分析, 使用了GIS地理数据库管理和空间对象叠置分析功能, 说明了在GIS环境中颗粒形状参数和属性参数的定义与获取方法; 对颗粒有限应变分析时, 使用了两种方法:最近邻弗莱法(TFry)和ASPAS法。以苏北盆地大丰—兴化勘探区块碎屑岩玻片为研究对象, 成功地分析了矢量化后的颗粒多边形形状及晶体塑性变形, 为构造地质分析和应用领域提供了一个方便快捷的工具。
关键词地理信息系统    形状参数    最近邻弗莱法    ASPAS法    有限应变    
文章编号:1006-6616(2009)03-0218-08
ANALYSES ON GRAIN SHAPE AND LIMITED RESPONSE USING GIS
LI Bo , ZHANG Ting-shan , WANG Zhan-lei , DONG Yin-lei     
School of Resources and Environment, SWPU, Chengdu 610500, China
Abstract: By geographical database management and space objects overlay analysis of GIS, this article makes a use of GIS in the analyses on grain shape and limited response of rocks, bringing about an approach of defining and getting grain parameters of the shape and attributes in GIS environment. In the analysis on particle limited response, two methods were used: TFry and ASPAS. Taking the thin sections of clastic rock from Dafeng-Xinghua exploration blocks in north Jiangsu Basin as research objects, an analysis has been made on the polygon shape and crystal plastic property of grains after vectorization, providing a powerful method for structural geology analysis and application.
Key words: GIS    shape parameter    TFry    ASPAS    limited response    
0 引言

岩石的变形程度通常以岩石的应变来度量。岩石有限应变测量就是利用岩石中某些标志体(如化石、鲕粒、矿物等)的形态、分布和物性来确定岩石的应变状态。岩石有限应变测量有助于认识岩石中构造的形成机制, 能够恢复岩石变形前的状态, 而且区域性的系统应变测量还可以建立区域的应变场。岩石有限应变测量不仅可以帮助我们认识岩石中构造的性质, 有助于查明这些构造的形成机制; 而且也是对岩石变形进行定量研究的基础, 是现代构造地质研究中不可缺少的组成部分[1]

对岩石有限应变的测量, 往往要收集大量构造要素数据, 对数据进行加工和整理, 然后绘制反映地质体构造几何学和运动学的各种构造地质图件, 并进行受力变形的点、面、体的应力和应变分析, 用手工操作要花费大量的人力和时间, 效率低下[2]。地理信息系统(GIS)是在计算机软硬件支持下, 对地理空间信息进行采集、存储、管理、运算、分析、显示和描述的技术系统, 非常适合于地球科学领域相关问题的研究, 为基础地质工作注入了新的生命力, 地理信息系统的应用使地质勘查工作出现了质的飞跃[3]。GIS的海量数据处理能力和使用形状分析得到颗粒形状信息的方便性, 给地理信息系统的应用带来了极大的优越性[4]。GIS在构造地质学中的应用, 可以简单地归纳为:数据信息存贮与分析; GIS中构造地质制图; GIS构造模拟与分析。

本文将GIS应用于岩石颗粒形状和有限应变分析, 为构造地质分析提供了一个方便快捷的工具。

1 形状参数与属性参数

本文主要对岩石颗粒形状的3个几何特征进行描述:粒度、形状和颗粒方向。一个颗粒的粒度由该类颗粒玻片下几何直径参数d来决定, 而某一样品的粒度取决于样品的2个参数:不同直径范围百分比面积和直径出现的频率。

颗粒的形状也使用2个参数进行描述, 即颗粒椭圆率R和形状参数S(颗粒面积周长平方)。颗粒椭圆率被定义为颗粒的最大内接椭圆的长轴长度a和短轴长度b的比值(a/b), 由于椭圆率对颗粒的整体形状敏感, 而对颗粒边界的细节形状不敏感, 所以椭圆率参数R是作为颗粒形状参数S的补充。S的计算公式为:

其中Ai是样品中单个颗粒面积, Pi是单个颗粒的周长。

从公式中可以看出:S与粒度无关, S=1时表明颗粒形状是圆形, 当颗粒形状与圆形偏差增大时S值减小。颗粒形状参数S到完美椭圆曲线E的距离可作为测量颗粒边界不规则程度的依据[5]

颗粒方向被定义为最大内接椭圆的长轴和基准线之间的夹角ϕ, 基准线方向可以根据研究的需要指定。角度ϕ和椭圆率RRf/ϕ应变分析方法的参数[6]

在ArcGIS软件中, 每个岩石颗粒作为多边形对象存储在空间数据库中, 每个多边形对象都有自己对应的形状接口(shape), 通过接口可以获得多边形的面积、周长和质心; 而多边形最大内接椭圆参数可以通过ArcToolbox中的Zonal Geometry函数生成(见表 1)。

表 1 多边形形状参数获得方式 Table 1 Methods of obtaining shape parameters of polygon

对于获取的颗粒属性数据, 可直接在GIS多边形对象属性表中建立子类(subtype)来进行分类, 把矢量化后的颗粒组分属性值分为石英、长石和岩屑等等, 便于以后对玻片命名时统计各组分含量。

2 GIS中岩石有限应变分析原理 2.1 最近邻弗莱法(TFry)

在许多地质情况下, 岩石颗粒点的初始分布不是随机的, 而是呈集合体散乱分布的, 但其间距离在统计上近于常数。这样分布的形成一般是因为物体具有特征的初始大小, 这就意味着最近邻的心对心距离是受堆积的几何约束所控制的[7]。这种集合体在沉积过程中变形后, 集合体中心对心间距是连线方向上伸长度的函数, 这就为测量应变提供了可能性。

弗莱法(Fry)基本思想:设有一统计上粒度均匀的颗粒(简化为二维等半径圆)集合体产生一套近邻中心, 对于中心点所有周围的点都具有反映平均粒径和堆积类型的空间分布。标志体中心分布如同无数半径不同的共心圆, 变形后这一共心圆变为共心椭圆, 这些椭圆就是在外力作用下的和应变椭圆具有相同形状的椭圆(见图 1)。

图 1 Fry原理示意图 Figure 1 Theory diagram of Fry method

由于GIS中矢量对象之间位置的拓扑关系也被保存在空间数据库中, 对象与对象的空间位置很容易获得, 所以在GIS中对均匀应力的应变椭圆求解时, 对Fry法稍加修改, 简化为只求邻接颗粒的中心点应变后所在位置, 即为TFry方法。

2.2 ASPAS方法

ASPAS方法示意图见图 2。该方法基于直线截断法理论[8], 将颗粒边界图放置在以参考点(P0)为中心、间隔角为Δϕ的放射线网络(ASPAS)上进行叠加。在叠加处理中, ASPAS中每条射线被颗粒边界截断分为数目不定且长度不等的线段。ASPAS中每条射线对颗粒多边形平均截断长度(d)是射线方向(ϕ)的函数, 因此, 我们可以画出表现函数d (ϕ)的玫瑰图。玫瑰图绘制取决于参考点Pi的位置, 在不同参考点(P0P5)位置上绘制的玫瑰图可以用来检验颗粒应变的局部效应。如果综合不同参考点绘制的玫瑰图来计算所有玫瑰图截断长度的加权平均, 可以生成完全ASPAS玫瑰图(total ASPAS graphy), 该图中放射线结尾点使用最小二乘法可以拟合出一个椭圆, 称这个椭圆为完全ASPAS椭圆, 在应变同质、均匀的情况下, 完全ASPAS椭圆在形状和方向上接近应变椭圆[9]

图 2 ASPAS方法示意图 Figure 2 Diagram of ASPAS method
2.3 应变椭圆参数求解(Rα)

TFry法和ASPAS方法在处理过程的最后都需要使用最小二乘法处理得到的点位图, 拟合出一个椭圆, 该椭圆可以近似看作是样品颗粒的应变椭圆, 其表达式为:

(1)

想要获得椭圆的长短轴和方向, 使用x轴的顺时针方向到椭圆长半轴之间的夹角α进行轴旋转变换, 变换后等式(1)变为:

(2)
(3)

利用不变性质, 可以获得:

(4)
(5)

解这个方程组, 可以获得A′和C′, 得出椭圆的长半轴和短半轴:

(6)
(7)

最终可以得到应变椭圆的椭圆度R=a/b, 椭圆方向ϕ=arctan (B (A -C))/2。

3 GIS中进行样品分析

本文所分析的2块薄片样品取自苏北盆地大丰—兴化勘探区块的丰参1井和丰探9井的钻井岩心。研究区内钻井岩心描述和薄片鉴定结果表明:探区泰州组一段地层和岩性受成岩作用和沉积作用的共同影响, 其岩性包含长石岩屑砂岩、长石石英砂岩、岩屑石英砂岩及一些粉砂岩。

3.1 颗粒形状参数求取

在GIS中进行颗粒形状和应变分析之前, 首先要对显微镜下的岩石薄片进行矢量化。矢量化后的颗粒数目和边界形状是GIS中处理的数据来源, 决定了形状和应变分析的结果。GIS软件提供了方便快捷的矢量化工具(Raster to Polygon), 直接将栅格图像转换为矢量多边形格式, 当薄片图像模糊或者颗粒边界需要人工判断时, 可以先在其他辅助软件(PhotoShop或R2V)中进行图像处理和人工描绘后, 再转化为GIS矢量文件。图 3a即为颗粒薄片矢量化后对应的颗粒多边形边界图层。

图 3 薄片1在GIS中分析结果 Figure 3 Analyzed results of thin section 1 in GIS

在GIS中颗粒的形状参数和颗粒的属性信息存储在地理数据库(Geodatabase)中。地理数据库是为了更好管理和使用地理要素数据, 而按照一定的模型和规则组合起来的存储空间数据和属性数据的容器, 在该数据库中能够按表 1的方法计算颗粒的形状参数并进行统计分析, 得到如图 3b所示的颗粒方向直方图、颗粒椭圆度R和形状参数S统计图以及颗粒直径直方图。从统计图中可以看出, 该样品颗粒直径d主要位于0.0625 ~ 0.1250 mm之间, 粒级属于极细砂; 颗粒椭圆度R主要分布区间为1 ~ 3, 平均值为1.804, 形状属于次圆形, 其中曲线E为不同椭圆度完美椭圆的形状参数S, 图 3b-(2)中的点到曲线E的距离可作为测量颗粒边界不规则程度的依据; 颗粒方位ϕ没有明显的规律性, 分布比较平均, 其中个数最大频率出现在170°~ 180°区间。结合颗粒多边形的空间形状数据和组分属性数据, 将薄片1样品命名为细—极细粒岩屑长石砂岩。

薄片2样品的颗粒方位ϕ有很明显的集中分布性(见图 4), 方位角主要位于90°~ 140°之间, 颗粒的排列有明显的定向性, 这种方向性可能是由于水流对颗粒的冲刷作用产生的, 颗粒在流动介质中沉积时, 其延伸方向往往平行于或垂直于水流的方向。因此如果知道岩心地层的产状方位, 就可以在GIS中求解古水流方向。

图 4 薄片2颗粒边界图及方向直方图 Figure 4 Boundary map and histogram of grains in thin section 2

用于颗粒方位的一般统计方法是矢量求和的方法。设定测量数据总量为n, 将0°~ 180°分为k组, 假设i组的颗粒数为ni, 方位角为ϕi, i =1, 2, 3, …, k, 则n =n1 +n2 +… +nk, 令xi =ni cos2ϕi, yi =ni sin2ϕi, 则x =x1 +x2 +…+xk, y =y1 +y2 +…+yk, 设平均方位角为ϕ, 则tan2ϕ =y/x。用这个方法求得薄片2的颗粒平均方位ϕ=124.7°。

3.2 GIS中实现TFry分析

GIS中实现TFry法步骤:

① 根据颗粒多边形形状和空间位置, 使用IArea接口或者Zonal Geometry工具求出每个多边形对象中心, 存储在多边形对象属性表中的字段value (5)和value (6)中。

② 以每个多边形对象中心为参考点, 求出其每个邻接多边形中心, 并进行相应的坐标平移投点, 最后绘制出TFry中空椭圆环点图(见图 3c中的蓝色点)。查找多边形邻接中心并平移的部分代码如下:

Do Until (pFeature Is Nothing), 循环查找图层中每个多边形pFeature

    Set pGeo =pFeature. Shape

    Set pPolygon =pGeo

    Set pFilter=New SpatialFilter

    With pFilter, 空间查询, 关系类型为Touch (两个多边形有公共边)

        Set.Geometry =pPolygon

        . GeometryField = “SHAPE”

        . SpatialRel =esriSpatialRelTouches

    End With

    Set pFeatureCursor1 =pFeatureClass. Search (pFilter, False)

    Set pFeature1 =pFeatureCursor1. NextFeature

    Do Until (pFeature1 Is Nothing)

        Print # 1, pFeature1.Value (5)-pFeature.Value (5), pFeature1. Value (6)- pFeature.Value (6)

      End If

      Set pFeature1 =pFeatureCursor1. NextFeature

    Loop

    Set pFeature =pFeatureCursor. NextFeature

Loop

③ 使用最小二乘法拟合应变椭圆, 并进行椭圆参数的求解。最终结果如图 3c。薄片1颗粒有限应变椭圆R=1.66, α=170°。

3.3 GIS中实现ASPAS分析

叠置分析是地理信息系统中用来提取空间隐含信息的方法之一, 将有关主题层组成的各个数据层面进行叠置产生一个新的数据层面, 其结果综合了原来2个或多个层面要素所具有的属性, 同时叠置分析不仅生成了新的空间关系, 且数据层的属性联系起来产生了新的属性关系。从原理上来说, 叠置分析往往涉及到逻辑交、逻辑并、逻辑差等的运算, 在ArcGIS中实现ASPAS方法实际上就是线对象图层和多边形图层的交集操作[10]

交集操作(Intersect)得到2个图层的交集部分, 并且原图层的所有属性将同时在得到的新的图层上显示出来。输入图层为矢量化后的颗粒多边形图层和ASPAS放射线图层, 相邻射线之间Δϕ=3°, 执行ArcGIS中交集操作后(如图 3d), ASPAS射线图层中每条射线被多边形图层中的每个多边形边界截断为多个线段(图 3d中蓝色选中线段), 将每个方向上的截断线段长度进行统计后绘制出P0 ~ P5局部ASPAS玫瑰图(对应图 3d中的P0 ~ P5), 然后将局部玫瑰图综合生成P完全ASPAS玫瑰图(对应图 3dP)。从图 3d中可以看出, 薄片1在不同方位进行ASPAS后的截线段长度是不同的, 体现了颗粒有限应变的局部不同性。而完全ASPAS玫瑰图是整个薄片颗粒应变分析的综合值, 对其进行最小二乘法拟合椭圆后得到参数R=1.52, α=163°。

4 结论和问题

薄片图像矢量化过程是本文分析的重要步骤, 如果颗粒边界矢量化有误或位置不准确, 则直接影响到分析数据和结果。颗粒图片分辨率低和色彩分级多是影响矢量化的重要因素, 分辨率低造成图片模糊, 正交光下色彩分级多造成颗粒内部颜色不均, 使得边界自动识别困难。所以要在辅助软件(PhotoShop或R2V)中做图像处理, 如锐化、色彩分级和重采样等等。另外一个造成颗粒边界识别困难的原因是颗粒的成因和成分, 单晶石英和长石颗粒很容易识别, 而多晶石英或燧石, 颗粒组成成分消光不一致, 自动识别会误认为是多个颗粒, 要在单片光下人为识别。

岩石的晶体塑性变形有2种机制:粒间变形和粒内变形[11]。粒间滑动的发生可能早于粒内变形, 在这种情况下, 粒内变形虽未发生, 但是岩石却已变形, 所以单纯分析粒内变形并未认识岩石的全貌, 况且单矿物的古应力值与岩石的变形应力值差别也很大。当晶体塑性的几何形态因变形的原始颗粒重结晶为新的亚颗粒而改变, 原始的颗粒构造可能完全消失, 这种情况TFry和ASPAS方法将不能解释有关应变程度的信息。

GIS可以胜任岩石颗粒形状分析和晶体塑性有限应变测量的任务, 通过栅格矢量化多边形对象, 能够在地理数据库中进行颗粒空间数据和属性数据的存储和管理, 并提供了大量的空间叠置分析和二次开发编程工具来分析数据, 得出TFry和ASPAS结果, 为构造地质分析和应用领域增添了一个方便快捷的工具。

薄片1应变分析结果表明, 使用TFry方法和ASPAS方法得到的应变椭圆参数基本相近。由于TFry方法依据邻接多边形中心的求取, 所以TFry方法更适用于碎屑岩中没有重结晶或者重结晶比较少的情况, 要求颗粒之间的接触情况为线接触; 而ASPAS方法则可以应用于早期重结晶碎屑岩颗粒中。

参考文献
[1]
郑亚东, 常志忠. 岩石有限应变测量及韧性剪切带[M]. 北京: 地质出版社, 1985, 1-185.
ZHENG Ya-dong, CHANG Zhi-zhong. Finite strain measurement and ductile shear zone[M]. Beijing: Geological Publishing House, 1985, 1-185. (in Chinese)
[2]
朱志澄. 构造地质学[M]. 武汉: 中国地质大学出版社, 2003.
ZHU Zhi-cheng. Structural geology[M]. Wuhan: China University of Geosciences Press, 2003. (in Chinese)
[3]
李跃辉. 地理信息系统(MapGis)在地质制图中的应用[J]. 地质力学学报, 2006, 12(2): 274-278.
LI Yue-hui. Application of geographical information system (MapGis) in geological map production[J]. Journal of Geomechanics, 2006, 12(2): 274-278. (in Chinese)
[4]
Bonham-Carter G F. Geographic information systems for geoscientists:Modelling with GIS[M]. Oxford: Pergamon Press, 1994.
[5]
Fernández F J, Mené ndez-Duarte R, Aller J, et al. Application of Geographical Information Systems to shape-fabric analysis[J]. Geological Society of London, Special Publications, 2005, 245: 409-420. DOI:10.1144/GSL.SP.2005.245.01.20
[6]
兰姆塞J G, 胡伯M I. 现代构造地质学方法, 第一卷, 应变分析[M]. 刘瑞, 常志忠, 张荣昌, 译. 北京: 地质出版社, 1991.
Ramsay J G, Huber M I. The techniques of modern structural geology, Volume 1, strain analysis[M]. LIU Rui-xun, CHANG Zhi-zhong, ZHANG Rong-chang, (Translated). Beijing: Geological Publishing House, 1991.
[7]
石磊, 黄维平, 孟宪刚. 岩石变形图像分析的计算机统计方法初探[J]. 地质力学学报, 1998, 4(3): 91-96.
SHI Lei, HUANG Wei-ping, MENG Xian-gang. Determination of rock deformation based on analysis of images of microstructures[J]. Journal of Geomechanics, 1998, 4(3): 91-96. (in Chinese)
[8]
Exner H E. Analysis of grain-and particle-size distribution in metallic materials[J]. International M etallurgical Reviews, 1972, 17(1): 25-42.
[9]
Launeau P, Robin P Y F. Fabric analysis using the intercept method[J]. Tectonophysics, 1996, 267(1-4): 91-119. DOI:10.1016/S0040-1951(96)00091-1
[10]
汤国安, 杨昕. ArcGIS地理信息系统空间分析实验教程[M]. 北京: 科学出版社, 2006, 1-480.
TANG Guo-an, YANG Xin. Spatial analysis by ArcGIS Geographic Information System[M]. Beijing: Science Press, 2006, 1-480. (in Chinese)
[11]
刘瑞. 显微构造地质学[M]. 北京: 北京大学出版社, 1988, 1-235.
LIU Rui-xun. Microstructural Geology[M]. Beijing: Peking University Press, 1988, 1-235. (in Chinese)