一种侵入岩体发育时序获取方法与流程

本发明涉及地理信息和地质学领域，尤其涉及一种侵入岩体发育时序获取方法。

背景技术：

岩体切割律又称穿插关系，就侵入岩与围岩的关系来说，总是侵入者年代新，被侵入者年代老，这就是切割律。这一原理还可被用来确定具交切关系或包裹关系的任何两地质体或地质界面的新老关系，即切割者新，被切割者老；包裹者新，被包裹者老。如侵入岩中捕虏体的形成年代比侵入体的老；砾岩中的砾石本身形成年代比砾岩的老；被断层切割的地层或火成岩形成的年代比断层的年代老。目前，岩体的发育时间先后顺序主要依靠专家利用切割律原理目视判别，效率低，且判别质量因人而异。

技术实现要素：

发明目的：本发明针对现有技术存在的问题，提供一种侵入岩体发育时序获取方法，可以有效提高判别质量和判别效率。

技术方案：本发明所述的侵入岩体发育时序获取方法包括：

(1)根据侵入岩体的地质剖面矢量面图层，构建岩体集合和岩体邻接矩阵；

(2)根据岩体邻接矩阵从岩体集合中读取任意两个邻接岩体sa、sb，并分别获取这两个岩体的分块集合pa和pb；

(3)基于两侧分布规则或半包围规则进行邻接岩体切割关系判别；

(4)循环执行步骤(2)-(3)，直至完成所有邻接岩体切割关系判别；

(5)根据岩体的切割关系，生成岩体切割关系矩阵；

(6)基于岩体切割关系矩阵，运用归并排序方法，生成岩体的发育时序。

进一步的，步骤(1)包括：

(1-1)加载侵入岩体的地质剖面矢量面图层，得到所有岩体集合s＝{sn|n＝1,2,…,sn}，sn表示第n个岩体，sn表示岩体数量；

(1-2)创建大小为sn*sn的岩体邻接矩阵；

(1-3)根据地质剖面矢量面图层，判定不同岩体间的邻接关系，当两个岩体邻接时，将对应的岩体邻接矩阵元素值赋值为1；否则，赋值为0。

进一步的，步骤(2)包括：

(2-1)按照从左向右、从上向下顺序，从岩体邻接矩阵中读取两个尚未读取过的岩体sa、sb；其中a、b表示岩体序号，且a、b∈{1,2,…,sn}；

(2-2)若sa、sb对应的岩体邻接矩阵元素值为1，则执行步骤(2-3)；否则返回执行步骤(2-1)；

(2-3)分别读取岩体sa、sb中所有的分块，形成对应的分块集合pa＝{aα|α＝1,2,…,an}和pb＝{bβ|β＝1,2,…,bn}，其中，aα表示岩体sa第α个分块，bβ表示岩体sb第β个分块，an、bn分别为sa、sb中所有的分块数量。

进一步的，步骤(3)包括：

(3-1)从分块集合pa中提取出邻接个数大于等于2的分块，存入子集sa＝{ai|i＝1,2,…,ai}，其中，ai表示子集sa中第i个岩体分块，ai为分块数量；

(3-2)从子集sa中读取任意一个元素ai；

(3-3)从分块集合pb中获取与ai邻接的分块，存入子集sb＝{bj|j＝1,2,…,bj}，其中，bj表示子集sb中第j个岩体分块，bj为邻接ai的分块个数；

(3-4)判断子集sb中任意两个分块bj、bj+1与ai是否满足两侧分布规则，若是，则判定sa切割sb，并执行步骤(3-7)，否则执行步骤(3-5)；

(3-5)归并分块集合pb，构建sb的整块岩体rb；

(3-6)判断ai与rb是否满足半包围规则，若是，则判定sa切割sb，并执行步骤(3-7)，否则，判定sa、sb切割关系未知，并执行步骤(3-7)；

(3-7)返回执行步骤(3-2)，直至子集sa中所有元素都被遍历。

进一步的，步骤(3-4)中两侧分布规则的判断方法的包括：

(3-4-1)获取bj、bj+1的外接矩形的中心点获取ai的外接矩形的中心点和长边线段的端点pta(xta,yta)、pwa(xwa,ywa)；

(3-4-2)根据下式，基于端点pta，pwa计算长边的斜率k：

(3-4-3)根据下式，获取过中心点的直线pl的直线方程；

(3-4-4)根据下式，计算中心点与直线pl的位置关系指标r：

(3-4-5)若r≤0，则表示bj、bj+1位于ai两侧，满足两侧分布规则，判定sa切割sb；如果r＞0，表示bj、bj+1位于ai同侧，不满足两侧分布规则。

进一步的，步骤(3-6)中半包围规则的判断方法的包括：

(3-6-1)获取ai的外接矩形的中心点和rb的外接矩形frb；

(3-6-2)从外接矩形的四个角点中，获取外接矩形的最大纵坐标值yamax、最小纵坐标值yamin、最大横坐标值xamax和最小横坐标值xamin；

(3-6-3)从外接矩形frb的四个角点中，获取外接矩形frb的最大纵坐标值yrmax、最小纵坐标值yrmin、最大横坐标值xrmax和最小横坐标值xrmin；

(3-6-4)根据下式计算用于判断外接矩形是否位于外接矩形frb内部的指标t；

t＝(yamax-yrmax)(yamin-yrmin)(xamax-xrmax)(xamin-xrmin)

(3-6-5)若t>0，表示外接矩形位于外接矩形frb内部，则判定sa切割sb；若t≤0，表示外接矩形未位于外接矩形frb内部，则执行步骤(3-6-6)；

(3-6-6)根据下式计算用于判断中心点是否位于外接矩形frb内部的指标u；

(3-6-7)若u>0，则表示中心点位于外接矩形frb内部，则判定sa切割sb；若u≤0，表示中心点未位于外接矩形frb内部，则判定sa、sb切割关系未知。

进一步的，步骤(6)具体包括：

(6-1)基于岩体切割关系矩阵，采用归并排序方法，按照从老到新的顺序对岩体进行排序；其中，新老的判断准则为：对任意两个岩体，切割的岩体为新，被切割的岩体为老；

(6-2)按照岩体的排序序号，生成对应岩体的发育时序；

(6-3)将生成的发育时序添加到对应岩体的时序属性timeid内。

有益效果：本发明与现有技术相比，其显著优点是：本发明是一种侵入岩体发育时序获取方法，可以有效提高判别效率和判别质量，获取的发育时序对于复杂侵入岩体发育过程表达与模拟具有重要的研究与应用价值。

附图说明

图1是本实施例中采用的地质剖面数据；

图2是本发明提供的流程图；

图3是本发明采用的两侧分布规则示意图；

图4是本发明采用的半包围规则示意图；

图5是本实施例的岩体时序判别结果图(1代表最老，4代表最新)。

具体实施方式

下面对本发明技术方案作进一步详细的说明，本实施例的实验数据采用的是《普通地质学》(夏邦栋，1983)教材中的图6-4地质剖面(图1)；并且，为了聚焦侵入岩体这一研究对象，仅保留了侵入岩体及其围岩。该实验数据采用的投影坐标系为wgs84。下面结合附图，并通过描述一个具体的实施例，来进一步说明。

如图2所示，本实施例提供的侵入岩体发育时序获取方法包括：

(1)根据侵入岩体的地质剖面矢量面图层，构建岩体集合和岩体邻接矩阵。

该步骤具体包括：

(1-1)加载侵入岩体的地质剖面矢量面图层，得到所有岩体集合s＝{sn|n＝1,2,…,sn}，sn表示第n个岩体，sn表示岩体数量；在本实施例中，sn＝4；

(1-2)创建大小为sn*sn的岩体邻接矩阵；

(1-3)基于arcgisengineapi，根据地质剖面矢量面图层，判定不同岩体间的邻接关系，当两个岩体邻接时，将对应的岩体邻接矩阵元素值赋值为1；否则，赋值为0。在本实施例中，构建的岩体邻接矩阵如表1所示。

表1岩体邻接矩阵

(2)根据岩体邻接矩阵从岩体集合中读取任意两个邻接岩体sa、sb，并分别获取这两个岩体的分块集合pa和pb。

该步骤具体包括：

(2-1)按照从左向右、从上向下顺序，从岩体邻接矩阵中读取两个尚未读取过的岩体sa、sb；其中a、b表示岩体序号，且a、b∈{1,2,…,sn}；在本实施例中，取a＝1、b＝2示例；

(2-2)若sa、sb对应的岩体邻接矩阵元素值为1，则执行步骤(2-3)；否则返回执行步骤(2-1)；例如，参考表1，本实施例中，s1、s2的岩体邻接矩阵元素值为1，执行步骤(2-3)；

(2-3)分别读取岩体sa、sb中所有的分块，形成对应的分块集合pa＝{aα|α＝1,2,…,an}和pb＝{bβ|β＝1,2,…,bn}，其中，aα表示岩体sa第α个分块，bβ表示岩体sb第β个分块，an、bn分别为sa、sb中所有的分块数量。本实施例中，当a＝1、b＝2时，an＝1，bn＝3。

(3)基于两侧分布规则或半包围规则进行邻接岩体切割关系判别。

该步骤具体包括：

(3-1)从分块集合pa中提取出邻接个数大于等于2的分块，存入子集sa＝{ai|i＝1,2,…,ai}，其中，ai表示子集sa中第i个岩体分块，ai为分块数量。本实施例中，当a＝1、b＝2时，ai＝1。

(3-2)从子集sa中读取任意一个元素ai。

(3-3)从分块集合pb中获取与ai邻接的分块，存入子集sb＝{bj|j＝1,2,…,bj}，其中，bj表示子集sb中第j个岩体分块，bj为邻接ai的分块个数。本实施例中，当i＝1时，bj＝3。

(3-4)判断子集sb中任意两个分块bj、bj+1与ai是否满足两侧分布规则，若是，则判定sa切割sb，并执行步骤(3-7)，否则执行步骤(3-5)。

其中，如图3所示，两侧分布规则的判断方法包括：

(3-4-1)基于arcgisengineapi获取bj、bj+1的外接矩形的中心点获取ai的外接矩形的中心点和长边线段的端点pta(xta,yta)、pwa(xwa,ywa)；在本实施例中，当i＝1,j＝1时，xta＝311.21229,yta＝-110.33909，xwa＝311.21229,ywa＝-304.38990；

(3-4-2)根据下式，基于端点pta，pwa计算长边的斜率k：

在本实施例中，当i＝1,j＝1时，k为无限大；

(3-4-3)根据下式，获取过中心点的直线pl的直线方程；

(3-4-4)根据下式，计算中心点与直线pl的位置关系指标r：

本实施例中，当i＝1,j＝1，k为无限大时，r＝-5944.97592；

(3-4-5)若r≤0，则表示bj、bj+1位于ai两侧，满足两侧分布规则，判定sa切割sb；如果r＞0，表示bj、bj+1位于ai同侧，不满足两侧分布规则。在本实施例中，当i＝1,j＝1时，r＝-5944.97592<0，b1、b2位于a1两侧，则s1分割s2。

(3-5)归并分块集合pb，构建sb的整块岩体rb。

(3-6)判断ai与rb是否满足半包围规则，若是，则判定sa切割sb，并执行步骤(3-7)，否则，判定sa、sb切割关系未知，并执行步骤(3-7)。

其中，如图4所示，半包围规则的判断方法的包括：

(3-6-1)基于arcgisengineapi获取ai的外接矩形的中心点和rb的外接矩形frb；

(3-6-2)从外接矩形的四个角点中，获取外接矩形的最大纵坐标值yamax、最小纵坐标值yamin、最大横坐标值xamax和最小横坐标值xamin；

(3-6-3)从外接矩形frb的四个角点中，获取外接矩形frb的最大纵坐标值yrmax、最小纵坐标值yrmin、最大横坐标值xrmax和最小横坐标值xrmin；

(3-6-4)根据下式计算用于判断外接矩形是否位于外接矩形frb内部的指标t；

t＝(yamax-yrmax)(yamin-yrmin)(xamax-xrmax)(xamin-xrmin)

(3-6-5)若t>0，表示外接矩形位于外接矩形frb内部，则判定sa切割sb；若t≤0，表示外接矩形未位于外接矩形frb内部，则执行步骤(3-6-6)；

(3-6-6)根据下式计算用于判断中心点是否位于外接矩形frb内部的指标u；

(3-6-7)若u>0，则表示中心点位于外接矩形frb内部，则判定sa切割sb；若u≤0，表示中心点pai未位于外接矩形frb内部，则判定sa、sb切割关系未知。

(3-7)返回执行步骤(3-2)，直至子集sa中所有元素都被遍历。

(4)循环执行步骤(2)-(3)，直至完成所有邻接岩体切割关系判别。

(5)根据岩体的切割关系，生成岩体切割关系矩阵。

其中，具体矩阵生成方法为：若sa切割sb，则对应元素值[sa,sb]赋值为“<”；若sa被sb切割，[sa,sb]赋值为“>”；若无法判断sa、sb的切割关系，则[sa,sb]赋值为“<>”。例如，在本实施例中，s1与s2的关系为s1分割s2，所以赋值“<“到[s1,s2]，最终构建的岩体切割关系矩阵如表2所示。

表2岩体切割关系矩阵

(6)基于岩体切割关系矩阵，运用归并排序方法，生成岩体的发育时序。

该步骤具体包括：

(6-1)基于岩体切割关系矩阵，采用归并排序方法，按照从老到新的顺序对岩体进行排序；其中，新老的判断准则为：对任意两个岩体，切割的岩体为新，被切割的岩体为老；

(6-2)按照岩体的排序序号，生成对应岩体的发育时序；本实施例中，归并排序的结果为s4、s2、s3、s1，即s1、s2、s3、s4的对应排序序号为4、2、3、1，则可以将s1、s2、s3、s4的发育时序定义为4、2、3、1，当然也可以定义为d、b、c、a等等，按照序号生成即可。

(6-3)将生成的发育时序添加到对应岩体的时序属性timeid内。

在本实施例中，侵入岩体时序判别结果如图5所示。本发明实施例中基于arcgisengineapi提供部分gis操作，相关步骤也可以使用supermap、arcgisobject等软件的api进行相应gis操作。

以上所揭露的仅为本发明一种较佳实施例而已，不能以此来限定本发明之权利范围，因此依本发明权利要求所作的等同变化，仍属本发明所涵盖的范围。