基于ArcGIS的Fry法自动实现.pdf

基于A r c G I S的F r y法自动实现收稿日期:2 0 2 3 0 4 2 4;修订日期:2 0 2 3 0 6 0 1;编辑:王敏作者简介:倪永进(1 9 8 8),男,山东昌乐人,工程师,从事构造地质研究工作;Em a i l:n i y o n g j i n 1 11 2 6.c o m*通讯作者:常杰(1 9 7 0),女,山东烟台人,工程师,主要从事地质矿产勘查及研究工作;Em a i l:c h a n g j i e d k y s h a n d o n g.c n倪永进1,常杰1*,田鹏飞2,宋香锁1,梁吉坡1,仵康林1(1.山东省地质科学研究院,山东 济南 2 5 0 0 1 3;2.山东省第一地质矿产勘查院,山东 济南 2 5 0 1 1 0)摘要:F r y法是一种基于应变标志体相对位置的有限应变测量方法,在构造地质研究领域广泛应用。传统的纸上作业方式过程繁琐,耗时漫长;已有的计算程序因编程语言不再使用、可视化程度低或主观性较强也已较少使用。本文基于A r c G I S平台,使用P y t h o n语言编制了新的F r y法应变测量程序。在点要素类上人工标记应变标志体的中心点后,运行程序即可快速获得F r y图解;交互输入空心椭圆半长轴、半短轴和旋转角的可能范围、步长后,利用点密度法可计算获得最佳拟合椭圆,并输出对应数据。本文通过2个实例,证明该程序操作简便,结果客观。关键词:F r y法;有限应变测量;A r c G I S;p y t h o n;点密度法中图分类号:P 2 0 8;T P 3 1 7 文献标识码:A d o i:1 0.1 2 1 2 8/j.i s s n.1 6 7 2 6 9 7 9.2 0 2 3.0 9.0 0 4引文格式:倪永进,常杰,田鹏飞,等.基于A r c G I S的F r y法自动实现J.山东国土资源,2 0 2 3,3 9(9):2 3 2 7.N I Y o n g j i n,CHANG J i e,T I AN P e n g f e i,e t a l.A u t o m a t i c I m p l e m e n t a t i o n o f F r y M e t h o d B a s e d o n A r c G I SJ.S h a n-d o n g L a n d a n d R e s o u r c e s,2 0 2 3,3 9(9):2 3 2 7.0 引言岩石有限应变测量系利用岩石内应变标志体的形态和分布来确定岩石的应变状态,构成了现代构造地质学研究的重要手段。反映岩石变形前和变形后几何形态差异的实体称为应变标志物1,常用的如化石、鲕粒等。近半个世纪以来,大量应变测量方法被提出23,其中一类关注应变标志物单体或集合体变形前后的形态或方向变化,如反向轮法、R f法、莫尔圆图解法37,另一类关注应变标志物的相对位置关系,如心对心法23。F r y89将心对心法发展为一种简便而优雅的测量方法,该方法通过平移统计应变标志物中心点,形成一个空心椭圆来代表应变椭圆。传统的F r y法利用透明纸进行手工制图,操作过程简单但经多次重复,耗时长、效率低,且在确定最终应变椭圆时主观性强、误差大1 01 2。为解决这一问题,过去3 0年中多种F r y法自动实现方案被设计出来,其中基于P A S C A L1 31 4、F o r t r a n1 2,1 5等语言由于设计年代较早,输入、输出以点的坐标数据为主,对可视化支持较差;基于O b j e c t i v e-C1 6语言、C o r e l D R AW平台1 7的仅解决了F r y图的自动生成,仍需人工确定最终的应变椭圆。W a l d r o n&W a l l a c e提出了一种客观拟合应变椭圆的方法1 0,即点密度法,该方法计算空心椭圆及与其具有相同面积的椭圆环中的点密度,并最小化这一比值来确定最佳拟合椭圆。相较于C o r e l D R AW、A I、P h o t o s h o p等作图软件,A r c G I S作为一种地理信息系统,具有良好的坐标和空间分析能力。虽然F r y法本身并不依赖地理坐标,但A r c G I S平台便捷的坐标提取和空间处理功能为本方法的自动制图和点密度运算提供了优良的工具选择。选择该平台的另一个重要原因是,A r c G I S平台提供了强大的p y t h o n接口,其A r c p y功能包可以经由p y t h o n语言调用而实现大量快速的重复计算。本文利用p y t h o n语言编制程序,调用A r c G I S平台的空间处理功能,自动绘制F r y图;并通过人机32第3 9卷第9期 山 东 国 土 资 源 2 0 2 3年9月交互,选择合理的应变椭圆参数范围,采用点密度法求解最佳拟合的应变椭圆。1 F r y法回顾在阐述自动化实现方法前,有必要回顾F r y法的基本原理和手工作图实现方法89。F r y法不考虑应变标志体本身的形态,只衡量标志体中心点相对位置3,其基本思想是:理想状态下,岩石截面可以简化为由半径为r的圆紧密堆积构成的集合体(图1)。在该堆积状态下,任意圆的圆心周围都有6个距离为2r的点(B、C等)、6个点的距离为(D、E等)构成一系列同心圆,岩石均匀变形后,同心圆变为同心椭圆。天然状态下,虽然应变标志体半径差异较大且堆积方式多变,但由于其在空间上的分布是随机的,经点的叠加仍可构成近于完整的椭圆。a变形前紧密堆积的颗粒圆及同心圆;b变形后紧密堆积的颗粒椭圆及同心椭圆图1 F r y法原理图(改自文献2)该方法操作过程中可能出现误差的主要有2处,即标志点的选择和空心椭圆的确定。前者主要由于样品表面不平整或亚颗粒的存在所造成,但基于大量的统计,这种误差的影响是微乎其微的。而空心椭圆作为F r y法最后确定应变的直接依据,通常由人为划定,主观性较强。大量工作讨论过客观自动识别空心椭圆的方法1 01 2,但也指出了这种自动识别的关键挑战-伪影1 8,即一个与真实空心椭圆存在交角的、范围较小的应变标志点空缺区。2 原理与方法2.1 原理基于A r c G I S的自动F r y法实现主体分为2个步骤:(1)自 动 生 成F r y图。本 文 中,将 照 片 导 入A r c G I S平台后,创建原始点要素类,采用目视解译,人工绘制标志体中心点位置。之后以a r c p y工具包读取原始点要素类中各点的坐标,并按照F r y法,自动绘制F r y图。(2)人机交互拟合最佳椭圆。由于其尺寸的局限性和显著的角度偏差,伪影的存在通常不会给人工确定空心椭圆造成困扰,只会影响缺乏直观感觉的自动算法;与之相反,不同观察者对空白椭圆轴比和旋转角的判断却可能存在较大差异,这是自动算法可以避免的。因此本文采用交互方法,先由人工依据第(1)步自动生成的F r y图解给出应变椭圆可能的参数取值范围,再由程序在该范围内进行检验,确定最佳参数并生成椭圆。最佳椭圆的自动检验标准采点计数密度法1 0。该方法首先构造一个内部椭圆和与之紧邻的外部椭圆环,二者中心点、轴比、旋转角和面积均相同(图2)。当拟合良好时,内部椭圆对应空心椭圆,包含最低密度的点;外部椭圆环对应点密环带,包含最高密度的点。因此可以通过在参数空间内的系统检验,使外-内密度比达到最大值。2.2 具体实现(1)将照片导入a r c m a p平台,并创建点要素类,以点要素标记标志物的中心点。形成原始点集po r i=pn(xn,yn)|n=1,2,3。42第3 9卷第9期 山 东 国 土 资 源 2 0 2 3年9月图2 点密度法原理图(改自文献1 0)(2)将所有点移动v1=(x1,y1),此时p1点移动至(0,0),记录其余点移动后的坐标,生成新的点位置;针对po r i数列中所有点重复上述过程,即v2=(x2,y2),v3=(x3,y3);得到新的F r y点集Pf r y=pm(xm,ym)|m=1,2,3n(n-1)。依据Pf r y生成新的点要素类并自动成图。(3)经观察步骤(2)生成的图像,输入空白椭圆半长轴a、半短轴b、旋转角k可能的取值范围r a n-g e a、r a n g e b、r a n g e k。以半长轴(a)、半短轴(b),生成2个椭圆,并旋转角度(k):椭圆1:x2a2+y2b2=1椭圆2:x22a2+y22b2=1 xryr=c o s(k)-s i n(k)s i n(k)-c o s(k)xy 2个椭圆的特点为轴比相同、主轴方向相同,后者面积为前者的2倍。分别计算在椭圆1内的点数Ni n和在椭圆2内但不在椭圆1内的点数No u t,计算后者的占比。C=No u t/Ni n 在r a n g e a、r a n g e b、r a n g e k内组合取值,直至遍历整个空间。由于内部椭圆1和外部圆环的面积相等,因此面积项被约掉,二者范围内点的个数比即为密度比。选取c取最大值时内部椭圆(椭圆1)的参数,即为空白椭圆参数。值得说明的是,在包含最佳取值且排除伪影取值的前提下,取值范围的给定不影响最终结果,只会对计算速度产生较大影响。3 应用实例与结果选择2个实例作为应用实例进行实验。实例1照片(图3 a)来自文献1 0,为加拿大英属哥伦比亚地区上K a z a组砂岩,碳酸盐胶结,分选良好。在照片中采集到了1 4 2个点,并生成了2 0 0 2 2个点构成的F r y图。为寻找最合适的椭圆,搜索空间如表1,在a b k空间中共搜索了5 4 0 8个点,最终得到的椭圆a=1.7,b=0.9,c=1.1 5(图3 c),椭圆内包含1个点,椭圆外环包含点7 9个,c值(椭圆外环点数占比)为9 8.7 3%。实例2照片(图3 b)来自文献1 2,为瑞典G a d-d e d e地区砾岩,颗粒支撑,砾石长轴方向近水平。该照片中采集到了6 8个点,并生成了由4 5 5 6个点构成的F r y图。注意目视该图中空心椭圆近水平(0),但在左上右下方向(1 3 0 附近)存在一个伪影区。为了寻找最合适的椭圆,在a b k空间中搜索了1 3 1 0 4个点(表1),最终得到的椭圆a=1 0.2,b=1.8,c=0(图3 d),椭圆内包含2个点,椭圆外环包含点6 1个,c值为9 6.7 2%。a砂岩镜下照片及中心点;b砾岩露头照片及中心点;c砂岩的F r y图解和最佳拟合的空心椭圆;d砾岩的F r y图和最佳拟合的空心椭圆图3 岩石镜下实测照片及F r y图表1 搜索空间实例实例1实例2参数abkabk起始值1.30.70.4 560.8-0.2 6终止值3.72.31.4 91 350.2 6步长0.20.20.0 20.30.20.0 252第3 9卷第9期 地 质 与 矿 产 2 0 2 3年9月表1中a、b、k分别为拟合椭圆的半长轴、半短轴和旋转角度。由于F r y法对绝对长度不敏感,仅关注a与b的比值,因此长度单位可以是任意的。k采弧度值,以指向右为0,顺时针旋转为负,逆时针旋转为正。4 讨论与结论4.1 F r y法空心椭圆及其与应变的关系通常认为F r y图空心椭圆以外应当有一个点密环带,表征了距离应变标志点最近的质点在变形后的分布。但从本文2个实例对比来看,实例1中的点密环带相对清晰,c值可达到9 8.7 3%;实例2中点密环带表现并不明显,c值也低于实例1。该环带的显著性可能受2种因素制约,即应变均匀程度和颗粒接触关系。早期观点认为空心椭圆代表应变标志体的平均应变,点密环带代表全岩总应变,只有空心椭圆而无点密环带的F r y图表示标志体初始分布不均匀或变形不均匀1 5,1 9。另一种观点则认为该环带的显著性可能与颗粒的接触关系有关,即其在颗 粒 支 撑 环 境 下 的 显 著 性 要 高 于 基 质 支 撑 环境1 71 8,2 0。这一现象也引发了F r y法局限性的讨论。尽管该方法被用在多种岩性的应变测量中,甚至是用以评估矿化空间的分布2 1,但应当认识到,只有当符合如下假设时空心椭圆才可被视为应变椭圆:(1)应变标志点应当足够多。(2)应变标志点的空间分布应当是随机,或者说是反聚类的。(3)岩石内的应变应当是均匀且可被观察到的(反例如晶格位错等)。4.2 未来的发展方向本方法通过人机交互,实现了基于A r c G I S的F r y法自动制图和空心椭圆提取。但仍有部分值得改进之处,如应变标志点的自动提取和最佳椭圆的快速拟合。A r c G I S平台强大的空间处理能力可以快速的识别照片(栅格数据)中的颗粒(面),并提取应变标志体(面)的中心点(点)。但在实际实验中,照片的识别普遍被干扰,且干扰源不可控,即不同照片的干扰因素差别巨大。典型的,对于露头或手标本尺度的照片,干扰因素可以包括露头或标本表面不平整造成的照片折线、应变标志物与基质的颜色对比过小或边界不清晰;对于显微照片,干扰因素可以报考亚颗粒的存在、双晶或多晶的边界以及不均匀消光等。近几年在计算机视觉和大数据地质学领域,已有大量工作尝试解决图像识别问题2 22 5,但距离本工作的实际应用尚有效率和准确率问题需要进一步解决:在几十至百余个点的数量级上,人工修改点与新建点的速度相差无几。本方法计算过程主要耗时之处在于遍历可能的长轴短轴旋转角组合,由于反复调用a r c p y的点面生成方法、空间连接函数和读取要素字段值,造成该步骤耗时过长。在本文提供的方法中,该步骤相当于生成一个由a b k和拟合度构成的四维空间,在该空间内,前三者形成一个四维曲面,通过遍历该曲面上的节点,寻找到拟合度轴上的最大值点。而在当前的机器学习领域,该步骤似乎可以通过梯度下降解决,即在该曲面上随机选择一点,计算该点的梯度,并向梯度方向移动一定距离,该距离与梯度的大小呈正比,通过该过程可避免计算全部曲面节点,进而大幅度提高计算效率。致谢:本文成文过程中,山东省地质科学研究院地质研究所全体同事提供了支持。参考文献:1 L I S L E R J.S t r a i n a n a l y s i s f r o m p o i n t f a b r i c p a t t e r n s:A n o b-j e c t i v e v a r i a n t o f t h e F r y m e t h o dJ.J o u r n a l o f S t r u c t u r a l G e-o l o g y,2 0 1 0,3 2(7):9 7 5 9 8 1.2 R AM S AY J G.F o l d i n g a n d F r a c t u r i n g o f R o c k s M.N e w Y o r k:M c G r a w-H i l l B o o k C o,1 9 6 7:5 6 8.3 R AM S AY J G,H u b e r M I.T h e T e c h n i q u e s o f M o d e r n S t r u c-t u r a l G e o l o g y,v o l.1:S t r a i n A n a l y s i s M.A C A D EM I C P R E S S,1 9 8 3:1 2 0.4 周继彬,曾佐勋.岩石有限应变测量反向轮法的计算机C S D软件设计J.地球科学,2 0 0 1(1):1 0 5 1 0 9.5 马维峰,于在平,常宏.岩石有限应变测量R f/法的一种改进方法J.西北地质,1 9 9 9(4):1 8 2 1.6 李云安.岩石有限应变测量的莫尔圆图解网法J.地质科技情报,1 9 9 3(1):5 3 5 8.7 何绍勋.变形岩石的应变分析J.地质与勘探,1 9 8 2(2):2 8.8 F R Y N.R a n d o m p o i n t d i s t r i b u t i o n s a n d s t r a i n m e a s u r e m e n t i n r o c k sJ.T e c t o n o p h y s i c s,1 9 7 9,6 0(1/2):8 9 1 0 5.9 HANNA S S,F R Y N.A c o m p a r i s o n o f m e t h o d s o f s t r a i n d e-t e r m i n a t i o n i n r o c k s f r o m s o u t h w e s t D y f e d(P e m b r o k e s h i r e)a n d a d j a c e n t a r e a sJ.J o u r n a l o f S t r u c t u r a l G e o l o g y,1 9 7 9,1(2):1 5 5 1 6 2.62第3 9卷第9期 山 东 国 土 资 源 2 0 2 3年9月1 0 WA L D R ON J,WA L L A C E K D.O b j e c t i v e f i t t i n g o f e l l i p s e s i n t h e c e n t r e t o c e n t r e(F r y)m e t h o d o f s t r a i n a n a l y s i sJ.J o u r n a l o f S t r u c t u r a l G e o l o g y,2 0 0 7,2 9(9):1 4 3 0 1 4 4 4.1 1 V I N T A B,S R I VA S T AVA D C.R a p i d e x t r a c t i o n o f c e n t r a l v a c a n c y b y i m a g e a n a l y s i s o f F r y p l o t sJ.J o u r n a l o f S t r u c-t u r a l G e o l o g y,2 0 1 2,4 0(7):4 4 5 3.1 2 Y EHUA S,WE N J I AO X.A s t a t i s t i c a l e x a m i n a t i o n o f t h e F r y m e t h o d o f s t r a i n a n a l y s i sJ.J o u r n a l o f S t r u c t u r a l G e o l-o g y,2 0 1 1,3 3(5):1 0 0 0 1 0 0 9.1 3 E R S L R V E A.N o r m a l i z e d c e n t e r t o c e n t e r s t r a i n a n a l y s i s o f p a c k e d a g g r e g a t e sJ.J o u r n a l o f S t r u c t u r a l G e o l o g y,1 9 8 8,1 0(2):2 0 1 2 0 9.1 4 E R S L E V E A,G E H.L e a s t s q u a r e s c e n t e r t o c e n t e r a n d m e a n o b j e c t e l l i p s e f a b r i c a n a l y s i sJ.J o u r n a l o f S t r u c t u r a l G e o l o g y,1 9 9 0,1 2(8):1 0 4 7 1 0 5 9.1 5 王涛,王晓霞.F R Y应变测量法的计算机模拟制图程序及应用J.地球科学与环境学报,1 9 8 8(3):6 1 6 6.1 6 曾荣杨,张夏永,李建.运用O b j e c t i v e C实现F r y法有限应变测量J.四川有色金属,2 0 1 8(3):4 7 5 0.1 7 韩阳光,颜丹平,李政林.在C o r e l D RAW平台上进行F r y法有限应变测量的新技术J.现代 地质,2 0 1 5,2 9(3):4 9 45 0 0.1 8 G E N I E R F,E P A R D J L.T h e F r y m e t h o d a p p l i e d t o a n a u g e n o r t h o g n e i s s:P r o b l e m s a n d r e s u l t sJ.J o u r n a l o f S t r u c t u r a l G e o l o g y,2 0 0 7,2 9(2):2 0 9 2 2 4.1 9 C R E S P I J M.S o m e g u i d e l i n e s f o r t h e p r a c t i c a l a p p l i c a t i o n o f F r y s m e t h o d o f s t r a i n a n a l y s i sJ.J o u r n a l o f S t r u c t u r a l G e o l-o g y,1 9 8 6,8(7):7 9 9 8 0 8.2 0 S C O T T R.P a t e r s o n,A c o m p a r i s o n o f m e t h o d s u s e d i n m e a s-u r i n g f i n i t e s t r a i n s f r o m e l l i p s o i d a l o b j e c t s J.J o u r n a l o f S t r u c t u r a l G e o l o g y,1 9 8 3,5(6):6 1 1 6 1 8.2 1 V E A R N C OMB E J R,V E A R N C OMB E S.T h e s p a t i a l d i s t r i-b u t i o n o f m i n e r a l i z a t i o n:A p p l i c a t i o n s o f F r y a n a l y s i sJ.E c o-n o m i c G e o l o g y,1 9 9 9,9 4(4):4 7 5 4 8 6.2 2 C HOUDHUR Y K R,ME E R E P A,MU L C HR ON E K F.A u-t o m a t e d g r a i n b o u n d a r y d e t e c t i o n b y C A S R GJ.J o u r n a l o f S t r u c t u r a l G e o l o g y,2 0 0 6,2 8(3):3 6 3 3 7 5.2 3 B A R T O Z Z I M,B OY L E A P,P R I O R D J.A u t o m a t e d g r a i n b o u n d a r y d e t e c t i o n a n d c l a s s i f i c a t i o n i n o r i e n t a t i o n c o n t r a s t i m a g e sJ.J o u r n a l o f S t r u c t u r a l G e o l o g y,2 0 0 0,2 2(1 1):1 5 6 91 5 7 9.2 4 李勃,张廷山,王占磊,等.基于G I S的颗粒形状和应变分析J.地质力学学报,2 0 0 9,1 5(3):2 1 8 2 2 5.2 5 崔敏,蔡佳.构造弱变形区自动有限应变测量:以石英砂岩为例J.地质力学学报,2 0 1 4,2 0(2):1 5 9 1 6 4.A u t o m a t i c I m p l e m e n t a t i o n o f F r y M e t h o d B a s e d o n A r c G I SN I Y o n g j i n1,CHANG J i e1,T I AN P e n g f e i2,S ONG X i a n g s u o1,L I ANG J i p o1,WU K a n g l i n1(1.S h a n d o n g I n s t i t u t e o f G e o l o g i c a l S c i e n c e s,S h a n d o n g J i n a n 2 5 0 0 1 3,C h i n a;2.N o.1 E x p l o r a t i o n I n s t i-t u t e o f G e o l o g y a n d M i n e r a l R e s o u r c e s,S h a n d o n g J i n a n 2 5 0 1 1 0,C h i n a)A b s t r a c t:T h e F r y m e t h o d i s a f i n i t e s t r a i n m e a s u r e m e n t m e t h o d b a s e d o n t h e r e l a t i v e p o s i t i o n o f s t r a i n m a r k e r s.I t i s w i d e l y u s e d i n t h e f i e l d o f s t r u c t u r a l g e o l o g y r e s e a r c h.T h e t r a d i t i o n a l p a p e r-b a s e d o p e r a-t i o n m e t h o d i s c u m b e r s o m e a n d t i m e-c o n s u m i n g.E x i s t i n g c o m p u t i n g p r o g r a m s a r e l e s s c o mm o n l y u s e d d u e t o t h e d i s c o n t i n u a t i o n o f p r o g r a mm i n g l a n g u a g e s,l o w l e v e l o f v i s u a l i z a t i o n,o r s t r o n g s u b j e c t i v i t y.B a s e d o n t h e A r c G I S p l a t f o r m,b y u s i n g P y t h o n l a n g u a g e,a n e w f r y m e t h o d s t r a i n m e a s u r e m e n t p r o g r a m h a s b e e n d e v e l o p e d.A f t e r m a n u a l l y m a r k i n g t h e c e n t e r p o i n t o f t h e s t r a i n i n d i c a t o r o n t h e p o i n t f e a t u r e c l a s s,t h e f r y d i a g r a m c a n b e q u i c k l y o b t a i n e d b y r u n n i n g t h e p r o g r a m.A f t e r i n t e r a c t i v e l y i n p u t t i n g t h e p o s s i b l e r a n g e a n d s t e p s i z e o f t h e s e m i m a j o r a x i s,s e m i m i n o r a x i s,a n d r o t a t i o n a n g l e o f t h e h o l l o w e l-l i p s e,t h e p o i n t d e n s i t y m e t h o d c a n b e u s e d t o c a l c u l a t e t h e b e s t f i t t i n g e l l i p s e a n d o u t p u t t h e c o r r e s p o n d-i n g d a t a.I n t h i s p a p e r,t h r o u g h t w o e x a m p l e s,i t i s p r o v e d t h a t t h e p r o g r a m i s e a s y t o o p e r a t e a n d t h e r e-s u l t s a r e o b j e c t i v e.K e y w o r d s:F r y m e t h o d;f i n i t e s t r a i n m e a s u r e m e n t;A r c G I S;p y t h o n;p o i n t-d e n s i t y m e t h o d72第3 9卷第9期 地 质 与 矿 产 2 0 2 3年9月