基于立体几何学的岩体结构面产状取值方法研究

陈柯宇，赵建军，李清淼，赖琪毅，代欣然，邓杰

(地质灾害防治与地质环境保护国家重点实验室(成都理工大学)，成都 610059)

1 简介

青藏高原地区地质构造强烈、板块快速隆升，孕育了大量高陡岩质边坡，地质人员现场难以通过地质罗盘测量岩层产状，室内在地质图上求解产状工作繁琐且不直观[1-4]。随着复杂计算机视觉、无人机摄影和地质调查领域发展结合的深入，通过航空倾斜摄影可以获取厘米级分辨率的三维影像，为上述问题的解决提供了远程解决的新思路[5-8]。三维空间影像技术是能获取目标物体表部空间几何特征，实现远距离、无接触测量，通过室内数据叠加便可获取精度高、密度大的三维点云数据[9]。通过三维空间影像提取测点的经纬度及高程信息，继而计算高陡岩质边坡的结构面产状，整个过程相对安全、便捷[10]。目前计算产状常采用向量关系法和拟合平面法[11-14]。

本文研究对象为高陡岩质边坡坡表松动岩体的结构面产状，不涉及深部岩体。通过高分辨航空遥感正射影像获取研究区域的带有位置信息的三维模型，与现场踏勘对照，确认待测量的结构面的拟合点，引入计算隐伏活动断层倾向、倾角计算公式，通过Python语言将计算方法程序化并建立简洁直观的用户窗口。

2 计算原理

基于立体几何学的三点问题求解，相比传统的作图法具有公式化、程序化的特点。三点法求解结构面产状的充要条件是已知位于同一结构面上的三个点间的平距、高差及两两连线的方位角。结构面倾角的计算公式是已知三点之间平距和高差为变量的代数表达式。结构面走向计算公式是已知三点之间平距、高差及三点连线的方位角为变量的代数表达式。

2.1 两点间距离和方位角计算方法

借助航空斜摄影测量技术建立的三维图像可以远程快速获取单点的地理信息(经纬度、高程)。两点的平距及两点连线的方位角无法直接从三维图像上获取，本文利用球面三角学相关公式进行推导并求解。

将地球近似看作一个球体，求解球面上A、B两点的距离和方位角的示意图如图1所示。对于地球上两点A、B，假设A为基准点，B为测点，其经纬度信息分别为Aj、Aw、Bj、Bw，θa、θb、θc为基准点A、测点B、北极点N其中两点与地心O夹角。θd为面AON与面BON夹角。

图1 两点求解距离和方位角示意图

求解基准点A与测点B与地心O之间夹角θc可以利用球面余弦公式：

cos(θc)=cos(θa)×cos(θb)+sin(θa)×sin(θb)×cos(θc)=cos(90°-Bw)×cos(90°-Aw)+sin(90°-Bw)×sin(90°-Aw)×cos(Bj-Aj)

(1)

(2)

根据球面正弦公式：

(3)

转换得：

(4)

根据基准点A与测点B的位置关系，测点相对于基准点的方位角如下：

(1)B点相对于A位于第一象限，即Bj≥Aj且Bw≥Aw时：AB方位角为θA；

(2)B点相对于A位于第二象限，即Bj

(3)B点相对于A位于第三、四象限，即Bw

根据半正矢公式，球面上A、B两点距离LAB为：

(5)

2.2 结构面产状计算方法

岩层层面与任一假想水平面的交线称走向线，层面上与走向线垂直并沿斜面向下所引的直线叫倾斜线。倾斜线在水平面上的投影所指示的方向称岩层的倾向。层面上的倾斜线和它在水平面上投影的夹角称为倾角。

三点的高度关系和位置关系直接影响到结构面倾角及倾向的公式推导。以三点高程不等的位置关系为例，如图2中(a)、(b)所示，A、B、C分别为结构面ABC三点，A点高程最高，B点高程次之，C点高程最低，DN为真北方向。存在水平面A’CE，A’为A点在水平面上投影，E为AB与水平面交点，CE即为走向线。AC上存在D与B高程一致，在平面ABC作DG⊥DB，DG与CE交于F。可得DG⊥CE，DG(DF)即为倾斜线。B’、D’、G’为B、D、G在水平面上投影。令∠G’FG为α，α∈(0，0.5Π)，∠A’D’B’为β，β∈(0，Π)，∠CA’E为γ，γ∈(0，Π)，AC之间高差为ha，AB高差为hb，AB间平距c，AC间平距为b，BC间平距为a，B’D’长度为d，BF前进方向为φBF，AB前进方向为φAB，AC前进方向为φAC,BC前进方向为φBC。倾斜线DF与其在水平面A’CE投影夹角α即为倾角，DF’前进方向φBF即为倾向。

图2 三点求结构面产状示意图

(6)

(7)

(8)

(9)

如图2(a)所示位置关系下，结构面ABC倾向φDF的计算公式如下：

φDF=φAC+β-90°

(10)

由公式(6)、(7)、(8)、(9)、(10)可知, 结构面ABC的倾角α计算公式仅与点A、B、C高程及两点间水平距离有关,与点间的方位关系无关。而结构面ABC的倾向不仅与三点的高程及两点间水平距离有关，还与三点的位置关系、两点连线方位角有关。对于选取的三点不同位置关系下倾角计算公式总结如表1。

表1 结构面倾向计算公式

3 计算方法程序化

Python(Guido van Rossum，1991)是一种广泛使用的编程语言，可以构建计算程序所需的的函数式、指令式、结构化、面向对象的编程。Numpy是python语言的一个扩展程序库，提供大量的数学函数库。Tkinter是Tk GUI工具包的标准python接口，是python建立图形用户界面的标准工具包。它也有其他的流行的GUI库来代替，例如wxPython、PyQt。

借助以上三者，可以较好地实现算法的本地化并可视化窗口，有效降低后续操作人员的负担，同时视觉上更加直观且便于理解。设计流程如下：(1)搭建Python运行环境，并加载所需的扩展程序库；(2)基于Numpy的底层函数将计算公式转化为用户自定义的函数，通过输入函数建立人机交互路径；(3)通过Tkinter提供的底层函数建立图形用户界面。最终图形界面如图3所示。

图3 程序操作界面及Tkinter底层函数

设计完成的程序包括三个大的部分。窗口，用户显示屏上的某矩形区域，包含标准框架和控件。小工具，指图形用户界面中构成应用程序的任何构建块，包括容器、按钮和对话框。框架，是复杂布局的基本组织单位，是一个矩形区域，可以包含其他小部件。

程序各区域功能如下：A：输入区，按照提示依次输入三点的位置信息；B：确认按钮，完成前一项数据输入并进入下一项；C：提示区，对当前所需的操作进行提示；D：计算按钮，完成数据输入后，点击按钮开始计算；E，显示区：显示产状计算结果；F：备注区，备注提示信息或进行数据记录。

4 数据与测试

以藏东某沟谷高位岩质陡壁为例，如图4所示，为研究陡壁上部危岩落石的冲击影响，需要测明陡壁上部控制结构面产状，以确定危岩落石方量。由于后缘岩壁陡峭，人工测量仅能在岩壁下部及两侧开展，为测量现场人员难以到达的区域的岩体结构面产状，本工程借助航空倾斜摄影技术和产状计算程序进行调查。在大规模测量前，对程序计算产状与现场罗盘测量的产状进行对比，验证程序计算结果的准确性，测试结果如表2所示。室内测量的产状与现场罗盘测量的产状整体一致，倾向平均误差为3°，倾角平均误差为2°。现场调查产状测量的倾向和倾角中误差分别为4°和3°，现场容许的最大偶然误差值为中误差的3倍。室内测量的产状值在工程误差容许范围内，可以满足实际需求。

表2 现场实测产状与室内计算产状对比表

图4 测点无人机照片(#a)和航空倾斜摄影图(#b)对比图

上述结论表明借助三维影像和推导的倾角、倾向公式以及可以快速、准确的获取结构面的产状信息。同时不同的产状计算方法都可以通过Python语言和Tkinter库建立较好的图形用户界面。对于批量的数据也可以通过Python的循环语句，接入编写的计算函数进行批量处理。

5 结论

本文使用立体几何学知识及航空倾斜摄影技术总结出岩体表层结构面的快速取值方法。

(1) 通过三维立体影像确定结构面上三点经纬度、高程，利用封装的倾角、倾向计算程序可以快速、准确的计算结构面产状。

(2) 本文使用Python编程语言将倾向、倾角计算方法通过图形用户界面的方式提供给用户，让算法在实际工程应用中更加简洁直观。对于已有的解析几何中的向量关系和最小二乘法拟合平面的产状算法也可以通过本文流程实现窗口化。

(3) 限于现阶段航空测量精度和三维地形图只能提供岩体表面空间信息的问题，通过程序测量产状的方法现阶段仍无法完全取代传统的现场罗盘测量。

本文的创新点在于基于立体几何学和Python语言完成了产状计算程序开发的全流程。