CN103207419B - 隧道岩层产状三维测量方法 - Google Patents

Abstract

本发明公开了一种隧道岩层产状三维测量方法，该方法以Google Earth提供的海量DEM、影像等地理信息数据为数据源，利用Google Earth COM API实现从Google Earth上获取岩层分界点数据，进行岩层面拟合，计算岩层产状，并通过标记语言实现岩层面以及岩层产状要素的可视化显示。基于Google earth平台隧道岩层产状的三维测量方法相比用地质罗盘进行野外勘测，不仅省时省力，而且扩展了勘测范围，对自然条件恶劣、岩层跨度范围大的区域，亦可方便测出其岩层产状要素，是对传统测量方法的一个很好补充。

Description

隧道岩层产状三维测量方法

技术领域

本发明涉及工程地质勘察技术领域，具体涉及隧道岩层产状三维测量方法。

技术背景

岩层产状测量是工程地质勘察一项重要的基础技术工作，为隧道工程地质断面图的编制提供基础数据。传统岩层产状测量工作是在野外采用罗盘、坡度仪进行现场测量，存在工作量大、整体精度低、受野外条件限制等缺陷。随着三维数字地球空间信息技术的发展，逐步实现了基于计算机和网络平台的大区域范围的地形、地貌、地物信息快速获取，为岩层产状的三维空间测量技术奠定了基础。20世纪60年代，国外利用航片测量产状的理论开始出现并应用（Berger Z, Williams T H L, Anderson D W. Geologicstereo mapping ofgeologic structure with SPOT satellite data[J] .AAPGBulletin, 1992,76, 76 (3) :101-102 .），但因受限于遥感影像的分辨率和处理能力而没有得到推广。20世纪90年代，开始技术人员尝试利用航空立体像确定岩层产状（陈建平,苗放. 利用航空立体像对确定岩层产状的计算机方法[J]. 国土资源遥感,1991,10(4):40-45.CHEN Jianping, MIAO Fang．The calculation method for layer attitude with airphoto pair[J]. Remote Sensing for Land & Resources, 1991, 10(4):40-45.），此类方法可以获得较高精度，但存在航空立体像的获取成本高、处理周期长等缺点，不适合进行大范围、批量产状计算。

发明内容

本发明的目的是针对上述技术问题，提供一种精度高、处理周期短的隧道岩层产状三维测量方法，同时该方法还适合进行大范围、批量产状计算。

为实现此目的，本发明所设计的隧道岩层产状三维测量方法，其特征在于，它包括如下步骤：

步骤1：在Google Earth（谷歌地球）软件数据库中人工选取多个待测地表岩层分界点的大地坐标，选取的待测地表岩层分界点的个数不少于3个；

步骤2：采用高斯投影将上述待测地表岩层分界点大地坐标转换为平面坐标；

步骤3：在步骤2的平面坐标系中任意选取三个点确定基准平面方程：

A₀X+B₀Y+C₀Z+1=0                  （1）

其中，A₀、B₀、C₀是基准平面方程系数，x，y，z是三维空间直角坐标；

步骤4：根据上述基准平面拟合出一个近似平面，即根据上述基准平面方程确定拟合平面方程，任意选取多个点所拟合的平面方程为AX+BY+CZ+1=0（2）

则拟合平面方程和基准平面方程之间满足如下关系：

A=A₀+A′，B=B₀+B′，C=C₀+C′，其中A'、B'、C'为误差；

步骤5：确定待测地表岩层分界点(X_i,Y_i,Z_i)的拟合平面方程，将A=A₀+A′，B=B₀+B′，C=C₀+C′代入步骤4中确定的拟合平面方程AX+BY+CZ+1=0，得到待测地表岩层分界点(X_i,Y_i,Z_i)的拟合平面方程

V_i=A′+B′Y_i+C′Z_i+A₀X_i+B₀Y_i+C₀Z_i+1=A′X_i+B′Y_i+C′Z_i+M_i（3）

其中，M_i是将测量点i的坐标代入基准平面方程所得的误差，为常数；

步骤6：根据最小二乘法，用待测地表岩层分界点(X_i,Y_i,Z_i)拟合计算岩层面方程，应使：

$S=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(A^{\′}{X}_i+B^{\′}{Y}_i+C^{\′}{Z}_i+M_i)}^2---\left(4\right)$

最小，为求得最小值，S分别对A′,B′,C′求导，并令其等于0，得：

$\left\{\begin{array}{c}\frac{d\left(\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{V}_i^2\right)}{d{A}^{\′}}=0,\\ \frac{d\left(\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{V}_i^2\right)}{d{B}^{\′}}=0,\\ \frac{d\left(\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{V}_i^2\right)}{d{C}^{\′}}=0,\end{array}\right.---\left(5\right)$

将式（3）代入式（5）展开，得到线性方程组：

$\left\{\begin{array}{c}{A}^{\′}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{X}_i^2+B^{\′}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{X}_{i}Y+C^{\′}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{X}_i{Z}_i+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{M}_i{X}_i=0,\\ A^{\′}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{X}_{i}Y+B^{\′}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{Y}_i^2+C^{\′}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{Y}_i{Z}_i+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{M}_i{Y}_i=0,\\ A^{\′}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{X}_i{Z}_i+B^{\′}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{Y}_i{Z}_i+C^{\′}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{Z}_i^2+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{M}_i{Z}_i=0.\end{array}\right.---\left(6\right)$

解线性方程组（6），可求得A′,B′,C′，将A′,B′,C′代入A=A₀+A′，B=B₀+B′，C=C₀+C′，拟合出岩层面方程为：

(A₀+A′)X+(B₀+B′)Y+(C₀+C′)Z+1=0（7）

其中，Z为岩层边界点的高程；

步骤7：得到待测岩层在空间分布状态的产状要素描述信息，上述要素包括岩层的走向、倾向和倾角，其中，待测岩层的走向为岩层面与水平面的交线；

待测岩层的倾向为岩层面上最大倾斜线在水平面上投影所指的方向，计算方法是在上述岩层走向的基础上加或减90°；

待测岩层的倾角为岩层面与水平面的夹角；

步骤8：获取待测隧道中线上任意一点的大地坐标和方向，

步骤9：根据待测地表岩层面向待测隧道中线延伸相交的位置，确定待测隧道与待测地表岩层面的分界线，其中待测隧道中线起点和终点的三维空间坐标分别为(X₁,Y₁,Z₁)和(X₂,Y₂,Z₂)，则待测隧道中线方程为：

$\frac{X-X_1}{X_2-X_1}=\frac{Y{-Y}_1}{Y_2-Y_1}=\frac{Z-Z_1}{Z_2-Z_1}---\left(10\right)$

求解方程组：

$\left\{\begin{array}{c}\frac{X-X_1}{X_2-X_1}=\frac{Y-Y_1}{Y_2-Y_1}=\frac{Z-Z_1}{Z_2-Z_1}\\ (A_0+A^{\′})X+(B_0+B^{\′})Y+(C_0+C^{\′})Z+1=0\end{array}\right.---\left(11\right)$

方程组的解即为待测地表岩层面与待测隧道的中线在三维空间中的交点，依次计算出各交点到线路中线上的投影里程，即得到待测隧道岩层分界点里程；

步骤10：在空间坐标系下对式（7）的岩层面方程进行可视化操作，在该可视化操作中须为其拟定一个边界，取步骤1中待测岩层分界点平面坐标（x、y）的最大值和最小值点作为边界点，记为（MaxX、MaxY）、（MaxX、MinY）、（MinX、MaxY）和（MinX、MinY）,连接上述四个点（MaxX、MaxY）、（MaxX、MinY）、（MinX、MaxY）和（MinX、MinY）所组成的平面，即为待测岩层面边界；将上述四个边界点数据代入式（7）的岩层面方程，得到边界点的高程Z；考虑到Google Earth软件采用的是大地坐标系，因此还须进行高斯投影坐标反算，得到待测岩层边界点的大地坐标；

步骤11：编辑Google Earth软件中的标记语言（Keyhole MarkupLanguage，缩写为KML）文件，在标记语言文件的“多边形”标签中添加步骤10中获得的待测岩层边界点的大地坐标数据，并在标记语言文件的“描述”标签中添加步骤7中得到的待测岩层的走向、倾向和倾角信息，同时，进行标记语言文件的实时生成与加载，完成计算结果的实时可视化显示，这样就实现了隧道岩层产状三维测量。

所述步骤7中岩层的走向的计算方法为：AX+BY+CZ+1=0平面与平面Z=0的交线，岩层走向的计算公式为：

$\mathrm{\φ}=\arccos (-B/\sqrt{A^2+B^2}).---\left(8\right)$

所述步骤7中岩层倾角通过如下公式计算：

$\mathrm{\θ}=\arccos (C/\sqrt{A^2+B^2+C^2}).---\left(9\right)$

所述步骤11中标记语言文件的实时生成与加载通过C#编程实现。

本发明以Google Earth提供的海量DEM（Digital Elevation Model，数字高程模型)、影像等地理信息数据为数据源，利用Google Earth COM API（谷歌地球软件组件对象模型的应用程序编程接口）实现从Google Earth上获取岩层分界点数据，进行岩层面拟合，计算岩层产状，并通过标记语言实现岩层面以及岩层产状要素的可视化显示。基于Google earth平台隧道岩层产状的三维测量方法相比用地质罗盘进行野外勘测，不仅省时省力，而且扩展了勘测范围，对自然条件恶劣、岩层跨度范围大的区域，亦可方便测出其岩层产状要素，是对传统测量方法的一个很好补充。利用Google Earth作为三维地理信息系统平台进行岩层产状测量，省去了DEM数据、影像数据的获取与校准过程，提高了测量效率。同时通过Google Earth COM API实现了与Google Earth的交互，不仅能快速批量获取岩层分界点，而且可实时编辑调整有偏差的点，用最小二乘法拟合岩层面，保证了岩层面的拟合精度与准确性，并利用标记语言实现了计算结果的三维可视化表达。

附图说明

图1为本发明的详细流程图。

具体实施方式

以下结合附图和实施例对本发明作进一步的详细说明：

如图1所述的隧道岩层产状三维测量方法，其特征在于，它包括如下步骤：

步骤1：在Google Earth软件数据库中人工选取多个待测地表岩层分界点的大地坐标，选取的待测地表岩层分界点的个数不少于3个；

步骤2：采用高斯投影将上述待测地表岩层分界点大地坐标转换为平面坐标；

步骤3：在步骤2的平面坐标系中任意选取三个点确定基准平面方程：

A₀X+B₀Y+C₀Z+1=0                  （1）

其中，A₀、B₀、C₀是基准平面方程系数，x，y，z是三维空间直角坐标，基准平面方程系数A₀、B₀、C₀可由上述任意选取的三个点的坐标求出；

步骤4：根据上述基准平面拟合出一个近似平面，即根据上述基准平面方程确定拟合平面方程，任意选取多个点所拟合的平面方程为AX+BY+CZ+1=0（2）

则拟合平面方程和基准平面方程之间满足如下关系：

A=A₀+A′，B=B₀+B′，C=C₀+C′，其中A'、B'、C'为误差；

步骤5：确定待测地表岩层分界点(X_i,Y_i,Z_i)的拟合平面方程，将A=A₀+A′，B=B₀+B′，C=C₀+C′代入步骤4中确定的拟合平面方程AX+BY+CZ+1=0，得到待测地表岩层分界点(X_i,Y_i,Z_i)的拟合平面方程V_i=A′X_i+B′Y_i+C′Z_i+A₀X_i+B₀Y_i+C₀Z_i+1=A′X_i+B′Y_i+C′Z_i+M_i（3）

其中，M_i是将测量点i的坐标代入基准平面方程所得的误差，为常数；

步骤6：根据最小二乘法，用待测地表岩层分界点(X_i,Y_i,Z_i)拟合计算岩层面方程，应使：

$S=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(A^{\′}{X}_i+B^{\′}{Y}_i+C^{\′}{Z}_i+M_i)}^2---\left(4\right)$

最小，为求得最小值，S分别对A′,B′,C′求导，并令其等于0，得：

$\left\{\begin{array}{c}\frac{d\left(\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{V}_i^2\right)}{d{A}^{\′}}=0,\\ \frac{d\left(\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{V}_i^2\right)}{d{B}^{\′}}=0,\\ \frac{d\left(\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{V}_i^2\right)}{d{C}^{\′}}=0,\end{array}\right.---\left(5\right)$

将式（3）代入式（5）展开，得到线性方程组：

$\left\{\begin{array}{c}{A}^{\′}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{X}_i^2+B^{\′}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{X}_{i}Y+C^{\′}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{X}_i{Z}_i+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{M}_i{X}_i=0,\\ A^{\′}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{X}_{i}Y+B^{\′}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{Y}_i^2+C^{\′}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{Y}_i{Z}_i+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{M}_i{Y}_i=0,\\ A^{\′}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{X}_i{Z}_i+B^{\′}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{Y}_i{Z}_i+C^{\′}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{Z}_i^2+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{M}_i{Z}_i=0.\end{array}\right.---\left(6\right)$

解线性方程组（6），可求得A′,B′,C′，将A′,B′,C′代入A=A₀+A′，B=B₀+B′，C=C₀+C′，拟合出岩层面方程为：

(A₀+A′)X+(B₀+B′)Y+(C₀+C′)Z+1=0             （7）

其中，Z为岩层边界点的高程；

步骤7：得到待测岩层在空间分布状态的产状要素描述信息，上述要素包括岩层的走向、倾向和倾角，其中，待测岩层的走向为岩层面与水平面的交线；

待测岩层的倾向为岩层面上最大倾斜线在水平面上投影所指的方向，计算方法是在上述岩层走向的基础上加或减90°；

待测岩层的倾角为岩层面与水平面的夹角；

步骤8：获取待测隧道中线上任意一点的大地坐标和方向，

步骤9：根据待测地表岩层面向待测隧道中线延伸相交的位置，确定待测隧道与待测地表岩层面的分界线，其中待测隧道中线起点和终点的三维空间坐标分别为(X₁,Y₁,Z₁)和(X₂,Y₂,Z₂)，则待测隧道中线方程为：

$\frac{X-X_1}{X_2-X_1}=\frac{Y{-Y}_1}{Y_2-Y_1}=\frac{Z-Z_1}{Z_2-Z_1}---\left(10\right)$

求解方程组：

$\left\{\begin{array}{c}\frac{X-X_1}{X_2-X_1}=\frac{Y-Y_1}{Y_2-Y_1}=\frac{Z-Z_1}{Z_2-Z_1}\\ (A_0+A^{\′})X+(B_0+B^{\′})Y+(C_0+C^{\′})Z+1=0\end{array}\right.---\left(11\right)$

方程组的解即为待测地表岩层面与待测隧道的中线在三维空间中的交点，依次计算出各交点到线路中线上的投影里程，即得到待测隧道岩层分界点里程；

步骤10：在空间坐标系下对式（7）的岩层面方程进行可视化操作，在该可视化操作中须为其拟定一个边界，取步骤1中待测岩层分界点平面坐标（x、y）的最大值和最小值点作为边界点，记为（MaxX、MaxY）、（MaxX、MinY）、（MinX、MaxY）和（MinX、MinY）,连接上述四个点（MaxX、MaxY）、（MaxX、MinY）、（MinX、MaxY）和（MinX、MinY）所组成的平面，即为待测岩层面边界；将上述四个边界点数据代入式（7）的岩层面方程，得到边界点的高程Z；考虑到Google Earth软件采用的是大地坐标系，因此还须进行高斯投影坐标反算，求出待测岩层边界点的经、纬坐标，从而得到待测岩层边界点的大地坐标；

步骤11：编辑Google Earth软件中的标记语言文件，在标记语言文件的“多边形”<LinearRing>标签中添加步骤10中获得的待测岩层边界点的大地坐标数据，并在标记语言文件的“描述”<description>标签中添加步骤7中得到的待测岩层的走向、倾向和倾角信息，同时，进行标记语言文件的实时生成与加载，完成计算结果的实时可视化显示，这样就实现了隧道岩层产状三维测量。

上述技术方案的步骤7中岩层的走向的计算方法为：AX+BY+CZ+1=0平面与平面Z=0的交线，岩层走向的计算公式为：

$\mathrm{\&phi;}=\arccos (-B/\sqrt{A^2+B^2}).---\left(8\right)$

上述技术方案的步骤7中岩层倾角通过如下公式计算：

$\mathrm{\&theta;}=\arccos (C/\sqrt{A^2+B^2+C^2}).---\left(9\right)$

上述技术方案的步骤11中标记语言文件的实时生成与加载通过C#编程实现。

下面以某隧道工程为实例，详细介绍以上各个步骤。

步骤1在Google earth软件中人工任意选取岩层分界线上的8个点，通过Google Earth COM API（谷歌地球软件组件对象模型的应用程序编程接口）获取其大地坐标；

步骤2：对上述岩层分界点，采用3度带高斯-克吕格投影，经过高斯投影变换，转换为高斯平面坐标，如表1所示。

表1 岩层分界点坐标

步骤3：采用最小二乘法式（6）对岩层分界点数据进行平面拟合计算，得到岩层面方程：

-3.2328×10^-7X+1.7728×10^-7Y-1.9171×10^-7Z+1=0        （13）

步骤4：根据式（8）、（9）和（13）计算出岩层产状：走向N61°E、倾角63°、倾向151°；

步骤5：计算边界点大地坐标，如表2;

表2 边界点的大地坐标

步骤6：在KML文件的<LinearRing>（多边形）标签中添加表2中的边界点大地坐标数据，在KML文件的<description>（描述）标签中添加步骤4中的岩层倾角、倾向、走向数据，最后Google earth软件中加载KML文件。

上述技术方案的步骤1具体来说是利用Google Earth提供的高分辨率卫星影像数据、地形数据，从宏观尺度上进行地表浅层岩层的判别，并确定岩层分界线。实现从Google Earth提取岩层分界点数据需要用到Google Earth COM API接口技术。其中IApplicationGE（谷歌地球应用程序接口）是Google Earth COM API的主要接口，它提供了从二维屏幕像素点坐标获取地形点坐标的函数，其函数原型为：GetPoint On Terrain From Screen Coords([in] double screen_x,[in] double screen_y,[out, retval] IPointOnTerrainGE**pPoint)，参数screen_x, screen_y为标准化的屏幕坐标，pPoint为函数返回的对应点的实际地理坐标。Google Earth地图窗口的左下角、右上角的坐标分别为（-1,-1），（1,1），窗口中心坐标为（0,0）。获取指定点坐标时，首先利用函数IApplieationGE::GetRenderHwnd获得Google Earth地图窗口句柄，然后利用这个句柄获得Google Earth地图窗口的长宽像素值，计算出具体点的标准化屏幕坐标。通过调用函数Get Point On Terrain From ScreenCoords([in] double screen_x,[in] double screen_y,[out, retval]IPoint OnTerrainGE **pPoint) 即可返回选取点pPoint的经度、纬度和高程值。在选点的过程中，如遇到选点不准或有偏移的情况，可通过重新拾取以及手动编辑的方式调整岩层分界点。

上述技术方案的步骤2中，采用高斯投影坐标正算法将待测地表岩层分界点大地坐标转换为平面坐标，具体公式为：

$\left\{\begin{array}{c}x=X+\frac{N}{2}\sin B\cos B{L}^2+\frac{N}{24}\sin B{\cos }^{3}B(5-t^2+9{\mathrm{\&eta;}}^2+4{\mathrm{\&eta;}}^4)L^4+\frac{N}{720}\sin B{\cos }^{5}B(61-58t^2+t^4)L^6\\ y=\mathrm{NconBL}+\frac{N}{6}{\cos }^{3}B(1-t^2+{\mathrm{\&eta;}}^2)L^3+\frac{N}{120}{\cos }^{5}B(5-18t^2+t^4+14{\mathrm{\&eta;}}^2-58t^2{\mathrm{\&eta;}}^2)L^5\end{array}\right.---\left(14\right)$

上述公式14为工程地质勘察领域常规算法公式，式中：B、L为椭球面上某点的大地坐标（弧度）;X、Y为某点投影后的平面直角坐标单位为m；X为纬度B对应的子午线弧长单位为m；N为所求点对应的卯酉圈子午线半径单位为m；t=tanB,η=cosB·e′，e′为第二偏心率。为避免y坐标出现负值，将x轴向西平移了500km，即所有点的y坐标值均加上500km。

上述技术方案的步骤10中，考虑到Google Earth软件采用的是大地坐标系，因此还须利用高斯投影变换坐标反算，求出岩层边界点的经、纬坐标（变换之前，y坐标须减去之前平移的500km）从而得到待测岩层边界点的大地坐标，利用经、纬坐标得到大地坐标是工程地质勘察领域常规方法。

本说明书未作详细描述的内容属于本领域专业技术人员公知的现有技术。

Claims (2)

1.一种隧道岩层产状三维测量方法，其特征在于，它包括如下步骤：

步骤1：在Google Earth软件数据库中人工选取多个待测地表岩层分界点的大地坐标，选取的待测地表岩层分界点的个数不少于3个；

步骤2：采用高斯投影将上述待测地表岩层分界点大地坐标转换为平面坐标；

步骤3：在步骤2的平面坐标系中任意选取三个点确定基准平面方程：

A₀X+B₀Y+C₀Z+1＝0   (1)

其中，A₀、B₀、C₀是基准平面方程系数，x，y，z是三维空间直角坐标；

步骤4：根据上述基准平面拟合出一个近似平面，即根据上述基准平面方程确定拟合平面方程，任意选取多个点所拟合的平面方程为AX+BY+CZ+1＝0   (2)

则拟合平面方程和基准平面方程之间满足如下关系：

A＝A₀+A′，B＝B₀+B′，C＝C₀+C′，其中A'、B'、C'为误差；

步骤5：确定待测地表岩层分界点(X_i,Y_i,Z_i)的拟合平面方程，将A＝A₀+A′，B＝B₀+B′，C＝C₀+C′代入步骤4中确定的拟合平面方程AX+BY+CZ+1＝0，得到待测地表岩层分界点(X_i,Y_i,Z_i)的拟合平面方程V_i＝A′X_i+B′Y_i+C′Z_i+A₀X_i+B₀Y_i+C₀Z_i+1＝A′X_i+B′Y_i+C′Z_i+M_i   (3)

其中，M_i是将测量点i的坐标代入基准平面方程所得的误差，为常数；

步骤6：根据最小二乘法，用待测地表岩层分界点(X_i,Y_i,Z_i)拟合计算岩层面方程，应使：

$S=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{(A^{\&prime;}{X}_i+B^{\&prime;}{Y}_i+C^{\&prime;}{Z}_i+M_i)}^2---\left(4\right)$

最小，为求得最小值，S分别对A′,B′,C′求导，并令其等于0，得：

$\left\{\begin{array}{c}\frac{d\left(\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{V_i}^2\right)}{{\mathrm{dA}}^{\&prime;}}=0,\\ \frac{d\left(\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{V_i}^2\right)}{{\mathrm{dB}}^{\&prime;}}=0,\\ \frac{d\left(\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{V_i}^2\right)}{{\mathrm{dC}}^{\&prime;}}=0.\end{array}\right.---\left(5\right)$

将式(3)代入式(5)展开，得到线性方程组： $\left\{\begin{array}{c}{A}^{\&prime;}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{X_i}^2+B^{\&prime;}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{X}_{i}Y+C^{\&prime;}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{X}_i{Z}_i+\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{M}_i{X}_i=0,\\ A^{\&prime;}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{X}_{i}Y+B^{\&prime;}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{Y_i}^2+C^{\&prime;}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{Y}_i{Z}_i+\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{M}_i{Y}_i=0,\\ A^{\&prime;}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{X}_i{Z}_i+B^{\&prime;}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{Y}_i{Z}_i+C^{\&prime;}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{Z_i}^2+\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{M}_i{Z}_i=0.\end{array}\right.---\left(6\right)$

解线性方程组(6)，可求得A′,B′,C′，将A′,B′,C′代入A＝A₀+A′，B＝B₀+B′，C＝C₀+C′，拟合出岩层面方程为：

(A₀+A′)X+(B₀+B′)Y+(C₀+C′)Z+1＝0   (7)

其中，Z为岩层边界点的高程；

步骤7：得到待测岩层在空间分布状态的产状要素描述信息，上述要素包括岩层的走向、倾向和倾角，其中，待测岩层的走向为岩层面与水平面的交线；

待测岩层的倾向为岩层面上最大倾斜线在水平面上投影所指的方向，计算方法是在上述岩层走向的基础上加或减90°；

待测岩层的倾角为岩层面与水平面的夹角；

步骤8：获取待测隧道中线上任意一点的大地坐标和方向，

步骤9：根据待测地表岩层面向待测隧道中线延伸相交的位置，确定待测隧道与待测地表岩层面的分界线，其中待测隧道中线起点和终点的三维空间坐标分别为(X₁,Y₁,Z₁)和(X₂,Y₂,Z₂)，则待测隧道中线方程为：

$\frac{X-X_1}{X_2-X_1}=\frac{Y-Y}{Y_2-Y_1}=\frac{Z-Z_1}{Z_2-Z_1}---\left(10\right)$

求解方程组：

$\left\{\begin{array}{c}\frac{X-X_1}{X_2-X_1}=\frac{Y-Y}{Y_2-Y_1}=\frac{Z-Z_1}{Z_2-Z_1}\\ (A_0+A^{\&prime;})X+(B_0+B^{\&prime;})Y+(C_0+C^{\&prime;})Z+1=0\end{array}\right.---\left(11\right)$

方程组的解即为待测地表岩层面与待测隧道的中线在三维空间中的交点，依次计算出各交点到线路中线上的投影里程，即得到待测隧道岩层分界点里程；

步骤10：在空间坐标系下对式(7)的岩层面方程进行可视化操作，在该可视化操作中须为其拟定一个边界，取步骤1中待测岩层分界点平面坐标(x、y)的最大值和最小值点作为边界点，记为(MaxX、MaxY)、(MaxX、MinY)、(MinX、MaxY)和(MinX、MinY),连接上述四个点(MaxX、MaxY)、(MaxX、MinY)、(MinX、MaxY)和(MinX、MinY)所组成的平面，即为待测岩层面边界；将上述四个边界点数据代入式(7)的岩层面方程，得到边界点的高程Z；考虑到GoogleEarth软件采用的是大地坐标系，因此还须进行高斯投影坐标反算，得到待测岩层边界点的大地坐标；

步骤11：编辑Google Earth软件中的标记语言文件，在标记语言文件的“多边形”标签中添加步骤10中获得的待测岩层边界点的大地坐标数据，并在标记语言文件的“描述”标签中添加步骤7中得到的待测岩层的走向、倾向和倾角信息，同时，进行标记语言文件的实时生成与加载，完成计算结果的实时可视化显示，这样就实现了隧道岩层产状三维测量。

所述步骤7中岩层的走向的计算方法为：AX+BY+CZ+1＝0平面与平面Z＝0的交线，岩层走向的计算公式为：

$\mathrm{\&phi;}=\arccos (-B/\sqrt{A^2+B^2}).---\left(8\right)$

所述步骤7中岩层倾角通过如下公式计算：

$\mathrm{\&theta;}=\arccos (C/\sqrt{A^2+B^2+C^2}).---\left(9\right)$

2.根据权利要求1所述的隧道岩层产状三维测量方法，其特征在于：所述步骤11中标记语言文件的实时生成与加载通过C#编程实现。