基于ArcGIS的Fry法自动实现

倪永进,常杰*,田鹏飞,宋香锁,梁吉坡,仵康林

(1.山东省地质科学研究院,山东 济南 250013;2.山东省第一地质矿产勘查院,山东 济南 250110)

0 引言

岩石有限应变测量系利用岩石内应变标志体的形态和分布来确定岩石的应变状态,构成了现代构造地质学研究的重要手段。反映岩石变形前和变形后几何形态差异的实体称为应变标志物[1],常用的如化石、鲕粒等。近半个世纪以来,大量应变测量方法被提出[2-3],其中一类关注应变标志物单体或集合体变形前后的形态或方向变化,如反向轮法、Rf-φ法、莫尔圆图解法[3-7],另一类关注应变标志物的相对位置关系,如心对心法[2-3]。Fry[8-9]将心对心法发展为一种简便而优雅的测量方法,该方法通过平移统计应变标志物中心点,形成一个空心椭圆来代表应变椭圆。

传统的Fry法利用透明纸进行手工制图,操作过程简单但经多次重复,耗时长、效率低,且在确定最终应变椭圆时主观性强、误差大[10-12]。为解决这一问题,过去30年中多种Fry法自动实现方案被设计出来,其中基于PASCAL[13-14]、Fortran[12,15]等语言由于设计年代较早,输入、输出以点的坐标数据为主,对可视化支持较差;基于Objective-C[16]语言、CorelDRAW平台[17]的仅解决了Fry图的自动生成,仍需人工确定最终的应变椭圆。Waldron &Wallace提出了一种客观拟合应变椭圆的方法[10],即点密度法,该方法计算空心椭圆及与其具有相同面积的椭圆环中的点密度,并最小化这一比值来确定最佳拟合椭圆。

相较于CorelDRAW、AI、Photoshop等作图软件,ArcGIS作为一种地理信息系统,具有良好的坐标和空间分析能力。虽然Fry法本身并不依赖地理坐标,但ArcGIS平台便捷的坐标提取和空间处理功能为本方法的自动制图和点密度运算提供了优良的工具选择。选择该平台的另一个重要原因是,ArcGIS平台提供了强大的python接口,其Arcpy功能包可以经由python语言调用而实现大量快速的重复计算。

本文利用python语言编制程序,调用ArcGIS平台的空间处理功能,自动绘制Fry图;并通过人机交互,选择合理的应变椭圆参数范围,采用点密度法求解最佳拟合的应变椭圆。

1 Fry法回顾

在阐述自动化实现方法前,有必要回顾Fry法的基本原理和手工作图实现方法[8-9]。Fry法不考虑应变标志体本身的形态,只衡量标志体中心点相对位置[3],其基本思想是:理想状态下,岩石截面可以简化为由半径为r的圆紧密堆积构成的集合体(图1)。在该堆积状态下,任意圆的圆心周围都有6个距离为2r的点(B、C等)、6个点的距离为(D、E等)……构成一系列同心圆,岩石均匀变形后,同心圆变为同心椭圆。天然状态下,虽然应变标志体半径差异较大且堆积方式多变,但由于其在空间上的分布是随机的,经点的叠加仍可构成近于完整的椭圆。

a—变形前紧密堆积的颗粒圆及同心圆;b—变形后紧密堆积的颗粒椭圆及同心椭圆

该方法操作过程中可能出现误差的主要有2处,即标志点的选择和空心椭圆的确定。前者主要由于样品表面不平整或亚颗粒的存在所造成,但基于大量的统计,这种误差的影响是微乎其微的。而空心椭圆作为Fry法最后确定应变的直接依据,通常由人为划定,主观性较强。大量工作讨论过客观自动识别空心椭圆的方法[10-12],但也指出了这种自动识别的关键挑战--伪影[18],即一个与真实空心椭圆存在交角的、范围较小的应变标志点空缺区。

2 原理与方法

2.1 原理

基于ArcGIS的自动Fry法实现主体分为2个步骤:

(1)自动生成Fry图。本文中,将照片导入ArcGIS平台后,创建原始点要素类,采用目视解译,人工绘制标志体中心点位置。之后以arcpy工具包读取原始点要素类中各点的坐标,并按照Fry法,自动绘制Fry图。

(2)人机交互拟合最佳椭圆。由于其尺寸的局限性和显著的角度偏差,伪影的存在通常不会给人工确定空心椭圆造成困扰,只会影响缺乏直观感觉的自动算法;与之相反,不同观察者对空白椭圆轴比和旋转角的判断却可能存在较大差异,这是自动算法可以避免的。因此本文采用交互方法,先由人工依据第(1)步自动生成的Fry图解给出应变椭圆可能的参数取值范围,再由程序在该范围内进行检验,确定最佳参数并生成椭圆。

最佳椭圆的自动检验标准采点计数密度法[10]。该方法首先构造一个内部椭圆和与之紧邻的外部椭圆环,二者中心点、轴比、旋转角和面积均相同(图2)。当拟合良好时,内部椭圆对应空心椭圆,包含最低密度的点;外部椭圆环对应点密环带,包含最高密度的点。因此可以通过在参数空间内的系统检验,使外-内密度比达到最大值。

图2 点密度法原理图(改自文献[10])

2.2 具体实现

(1)将照片导入arcmap平台,并创建点要素类,以点要素标记标志物的中心点。形成原始点集pori={pn(xn,yn)|n=1,2,3……}。

(3)经观察步骤(2)生成的图像,输入空白椭圆半长轴a、半短轴b、旋转角k可能的取值范围rangea、rangeb、rangek。以半长轴(a)、半短轴(b),生成2个椭圆,并旋转角度(k):

2个椭圆的特点为轴比相同、主轴方向相同,后者面积为前者的2倍。分别计算在椭圆1内的点数Nin和在椭圆2内但不在椭圆1内的点数Nout,计算后者的占比。

C=Nout/Nin

在rangea、rangeb、rangek内组合取值,直至遍历整个空间。由于内部椭圆1和外部圆环的面积相等,因此面积项被约掉,二者范围内点的个数比即为密度比。选取c取最大值时内部椭圆(椭圆1)的参数,即为空白椭圆参数。

值得说明的是,在包含最佳取值且排除伪影取值的前提下,取值范围的给定不影响最终结果,只会对计算速度产生较大影响。

3 应用实例与结果

选择2个实例作为应用实例进行实验。实例1照片(图3a)来自文献[10],为加拿大英属哥伦比亚地区上Kaza组砂岩,碳酸盐胶结,分选良好。在照片中采集到了142个点,并生成了20022个点构成的Fry图。为寻找最合适的椭圆,搜索空间如表1,在a-b-k空间中共搜索了5408个点,最终得到的椭圆a=1.7,b=0.9,c=1.15(图3c),椭圆内包含1个点,椭圆外环包含点79个,c值(椭圆外环点数占比)为98.73%。

表1 搜索空间

a—砂岩镜下照片及中心点;b—砾岩露头照片及中心点;c—砂岩的Fry图解和最佳拟合的空心椭圆;d—砾岩的Fry图和最佳拟合的空心椭圆

实例2照片(图3b)来自文献[12],为瑞典Gaddede地区砾岩,颗粒支撑,砾石长轴方向近水平。该照片中采集到了68个点,并生成了由4556个点构成的Fry图。注意目视该图中空心椭圆近水平(0°),但在左上—右下方向(130°附近)存在一个伪影区。为了寻找最合适的椭圆,在a-b-k空间中搜索了13104个点(表1),最终得到的椭圆a=10.2,b=1.8,c=0(图3d),椭圆内包含2个点,椭圆外环包含点61个,c值为96.72%。

表1中a、b、k分别为拟合椭圆的半长轴、半短轴和旋转角度。由于Fry法对绝对长度不敏感,仅关注a与b的比值,因此长度单位可以是任意的。k采弧度值,以指向右为0,顺时针旋转为负,逆时针旋转为正。

4 讨论与结论

4.1 Fry法空心椭圆及其与应变的关系

通常认为Fry图空心椭圆以外应当有一个点密环带,表征了距离应变标志点最近的质点在变形后的分布。但从本文2个实例对比来看,实例1中的点密环带相对清晰,c值可达到98.73%;实例2中点密环带表现并不明显,c值也低于实例1。该环带的显著性可能受2种因素制约,即应变均匀程度和颗粒接触关系。早期观点认为空心椭圆代表应变标志体的平均应变,点密环带代表全岩总应变,只有空心椭圆而无点密环带的Fry图表示标志体初始分布不均匀或变形不均匀[15,19]。另一种观点则认为该环带的显著性可能与颗粒的接触关系有关,即其在颗粒支撑环境下的显著性要高于基质支撑环境[17-18,20]。

这一现象也引发了Fry法局限性的讨论。尽管该方法被用在多种岩性的应变测量中,甚至是用以评估矿化空间的分布[21],但应当认识到,只有当符合如下假设时空心椭圆才可被视为应变椭圆:

(1)应变标志点应当足够多。

(2)应变标志点的空间分布应当是随机,或者说是反聚类的。

(3)岩石内的应变应当是均匀且可被观察到的(反例如晶格位错等)。

4.2 未来的发展方向

本方法通过人机交互,实现了基于ArcGIS的Fry法自动制图和空心椭圆提取。但仍有部分值得改进之处,如应变标志点的自动提取和最佳椭圆的快速拟合。

ArcGIS平台强大的空间处理能力可以快速的识别照片(栅格数据)中的颗粒(面),并提取应变标志体(面)的中心点(点)。但在实际实验中,照片的识别普遍被干扰,且干扰源不可控,即不同照片的干扰因素差别巨大。典型的,对于露头或手标本尺度的照片,干扰因素可以包括露头或标本表面不平整造成的照片折线、应变标志物与基质的颜色对比过小或边界不清晰;对于显微照片,干扰因素可以报考亚颗粒的存在、双晶或多晶的边界以及不均匀消光等。近几年在计算机视觉和大数据地质学领域,已有大量工作尝试解决图像识别问题[22-25],但距离本工作的实际应用尚有效率和准确率问题需要进一步解决:在几十至百余个点的数量级上,人工修改点与新建点的速度相差无几。

本方法计算过程主要耗时之处在于遍历可能的长轴—短轴—旋转角组合,由于反复调用arcpy的点—面生成方法、空间连接函数和读取要素字段值,造成该步骤耗时过长。在本文提供的方法中,该步骤相当于生成一个由a-b-k和拟合度构成的四维空间,在该空间内,前三者形成一个四维曲面,通过遍历该曲面上的节点,寻找到拟合度轴上的最大值点。而在当前的机器学习领域,该步骤似乎可以通过梯度下降解决,即在该曲面上随机选择一点,计算该点的梯度,并向梯度方向移动一定距离,该距离与梯度的大小呈正比,通过该过程可避免计算全部曲面节点,进而大幅度提高计算效率。

致谢：本文成文过程中,山东省地质科学研究院地质研究所全体同事提供了支持。