CN108253938A - Tbm破岩矿渣数字近景摄影测量识别及反演方法 - Google Patents

Abstract

一种TBM破岩矿渣数字近景摄影测量识别及反演方法，包括以下步骤：图像标定校正；TBM盾构机破岩产生矿渣；图像采集、姿态解算；图像处理、识别及分类；对已分类的矿渣进行室内物理力学实验得到矿渣对应的围岩岩性特征参数；建立图像信息和特征参数之间的矿渣关系型数据库；在矿渣关系型数据库搜索相似的图像信息，据此推测现场矿渣对应的围岩岩性特征参数；建立TBM掘进推力和掘进速度预测模型，预测前方TBM掘进推力和速度。本发明方法通过建立矿渣图像信息和矿渣对应的围岩岩性特征参数之间的矿渣关系型数据库，预测前方TBM掘进推力和速度，为推进系统设计及施工过程中参数实时调整提供重要依据，提高隧洞掘进开挖效率，降低施工成本。

Description

TBM破岩矿渣数字近景摄影测量识别及反演方法

技术领域

本发明涉及一种TBM破岩矿渣数字近景摄影测量识别及反演方法，具体涉及隧洞掘进机(TBM)、数字近景摄影测量技术和硬岩TBM设备掘进推力和速度的计算，尤其涉及TBM排渣系统、影像特征识别、分类、提取与特征匹配的方法，矿渣反射率、纹理、颜色信息数据采集及矿渣关系型数据库的建立以及TBM设备掘进推力和速度计算模型的建立。

背景技术

全断面隧道掘进机(TBM：Tunnel Boring Machine)广泛用于长隧洞机械化施工中。圆盘滚刀用于切削顽石和岩石，比较大的破碎岩石和顽石则由带状螺旋器排出装置运走,出渣一般采用连续胶带输送机或有轨运输矿车两种方式。而目前TBM破岩产生的矿渣仍然依靠现场人员经验来分辨，如何高效率地通过辨别矿渣来推测岩石的参数对前方TBM掘进推力和速度做出预测，从而为推进系统设计及施工过程中参数实时调整提供重要依据，一直是此领域技术人员探索的课题。

发明内容

本发明的目的在于提供一种TBM破岩矿渣数字近景摄影测量识别及反演方法，以此来匹配岩石的参数，对前方TBM掘进推力和速度做出预测，从而为推进系统设计及施工过程中参数实时调整提供重要依据。

为实现上述目的，本发明采取下述技术方案：

一种TBM破岩矿渣数字近景摄影测量识别及反演方法，其特征在于，包括以下步骤：

步骤1，图像标定校正：主要是由标定获得相机内参数矩阵和外参数矩阵，从而建立世界坐标系和摄像机坐标系以及图像坐标系的映射关系，再通过反变换以及灰度插值实现畸变图像的校正；

步骤2，TBM盾构机破岩产生矿渣；

步骤3，图像采集、姿态解算；

图像采集：在TBM破岩矿渣输送带的左边、右边和正上方的转台上分别固定摄像装置，所述摄像装置组成联机摄像系统，从输送带的左、右、上方采集矿渣图像；

姿态解算：采用PNP算法进行空间姿态计算，依据任意2D点的可视角度和对应3D点的角度一样，由一组2D点的映射中估计物体的3D姿态；

步骤4，图像处理、识别及分类；

图像处理：包括对矿渣图像进行灰度化处理，将彩色矿渣图像变成灰度矿渣图像；

图像识别：对灰度矿渣图像进行岩石图像灰度均值和纹理特征提取；计算多幅岩石图像的灰度均值，通过分析得到不同岩石图像灰度均值的范围，根据不同范围的图像灰度均值划分为不同种类的岩石；采用灰度共生矩阵进行纹理特征提取；

图像分类：采用不断训练的BP神经网络结构对矿渣进行分类；岩石灰度均值、纹理特征是作为BP神经网络结构的输入值数据，用来训练区分不同的岩石图像；图像分类的处理对象是岩石图像的灰度均值和纹理特征等提取的图像信息，处理结果是根据不同的图像信息分类矿渣图像，从而识别对应的矿渣；

步骤5，通过对矿渣图像的分类区分出所对应的矿渣种类，对已分类的矿渣进行室内物理力学实验，得到矿渣对应的围岩岩性特征参数，包括单轴抗压强度σ_c、围岩岩性类别R₁、泊松比μ、变形模量E_b、裂隙发育程度分级C_r；

步骤6，建立图像信息和特征参数之间的矿渣关系型数据库；所述图像信息是指矿渣图像的灰度均值和纹理特征的数字化信息，所述特征参数是指矿渣对应的围岩岩性特征参数；

步骤7，在矿渣关系型数据库搜索相似的图像信息，若存在相似图像信息，将矿渣参数上传至服务器主机储存的关系型数据库中，并据此推测现场矿渣对应的围岩岩性特征参数；若不存在相似图像信息，将新的图像信息加入矿渣关系型数据库中，更新矿渣关系型数据库；

步骤8，建立TBM掘进推力和掘进速度预测模型，预测前方TBM掘进推力和速度；TBM掘进总推力F预测模型如下：

式中，N为刀盘上装有的滚刀数量；π为圆周率；ω为滚刀刀圈宽度；r为滚刀半径；d为贯入深度；p₀为位于刀具正下方的基准压力；n＝0,1,2,…,∞；ψ为接触压力分布指数，取决于滚刀和岩石的参数；δ是滚刀岩石接触角；D是刀盘直径；l是护盾长度；是护盾与洞壁接触范围；p是接触压力，由护盾与地层接触情况决定；₁是护盾与岩石的摩擦系数，μ₂是后续设备与轨道间的摩擦系数；m是TBM后续设备质量，g是重力加速度；

TBM掘进速度预测模型如下：

Y(PR)＝β₀+β₁(σ_c)+β₂(R₁)+β₃(μ)+β₄(E_b)+β₅(C_r)

其中，Y(PR)为因变量，为掘进速度；β₀、β₁、β₂、β₃、β₄、β₅为各项自变量对应的系数；根据收集的现场试验数据，代入上述的正则方程，经过多元线性回归计算得到模型参数估计结果，得到的系数矩阵为：

将步骤7中推测得到的现场矿渣对应的围岩岩性特征参数单轴抗压强度σ_c、围岩岩性类别R₁、泊松比μ、变形模量E_b、裂隙发育程度分级C_r代入上述TBM掘进速度计算模型中即得到掘进速度Y(PR)。

上述步骤8中，将最新得到的R₁、σ_c、μ、E_b、C_r现场试验数据加入之前的数据表，从而不断通过多元线性回归修正并更新TBM掘进速度计算模型，并且将所述模型用于下一回合TBM掘进速度的预测，不断进行数据回归、预测速度、更新模型的循环。

上述步骤1中，图像坐标系和世界坐标系之间的转换公式如下：

其中，P＝K[R t]为相机矩阵或投影矩阵；[R t]为相机外参数矩阵，包括旋转矩阵R和平移向量t；L为相机内参数矩阵，x₀，y₀为主点坐标；f_x，f_y分别为图像x轴和y轴上的尺度因子，γ是描述两图像坐标轴倾斜程度的参数；Z_c为像空间坐标；x，y为像点的像平面坐标；X,Y,Z为物方点的物方空间坐标。

上述步骤1中，x₀，y₀，f_x，f_y计算公式如下：

式中：x₀，y₀为主点坐标；f_x，f_y分别为图像x轴和y轴上的尺度因子；l_i是相机内外方为元素的表达式，i＝1，2，3……11。

上述步骤1中，摄像装置的检校参数为：主距f像主点坐标(x₀,y₀)；畸变系数(k₁,k₂)为径向畸变系数；(P₁,P₂)为切向畸变系数；得到内参数矩阵L和畸变参数k_c；使用定标后的摄像装置内参数，求得的外方位元素，物方坐标真值依据共线条件方程推算像方点坐标；将实际量测的像方坐标经畸变改正后得到中误差像素大小。

上述步骤3的图像采集步骤，设置球铰调整摄像装置摄像角度，设置照明装置提供光照使所成图像清晰可辨，设置冲洗装置对摄像头进行清洗。

上述步骤3的姿态解算步骤，具体实现如下：

通过相机标定来获取相机内参数即基础矩阵和畸变向量，已知物体三维空间坐标系的点以及对应的图像像素坐标点；相机的参数矩阵L已通过相机标定得到，故可用PNP算法进行空间姿态计算，其依据任意2D点的可视角度和对应3D点的角度一样；采用P4P算法，取4个控制点：

一组2D和3D的对应点其中是单位方向矢量，通过标定矩阵L的逆将2D像素测量值变换为单位模长的3D方向矢量，得到：

未知量为从相机中心c到3D点p_i的距离d_i，其中由三角形Δ(c,p_i,p_j)余弦定理得：

式中：组合其它约束(f_ij，f_ik，f_jk)，并消去d_i,d_j可得如下关于的公式：

给定n个对应点，可以产生(n-1)(n-2)/2个三元组来得到的值的线性估计；可以用连续两个估计的比率来估计，然后通过平均来获得d_i；确定d_i后，可得到一系列尺度自由的点然后通过3D点之间的刚体变换，可将其变换到对应的样本点{p_i}，从而完成线性投影矩阵确定。

上述步骤4中的图像处理，在图像灰度化处理后，还进行了图像锐化和图像平滑处理：

图像灰度化，具体实现如下：

g(i,j)＝R＝G＝B＝0.2989R+0.5870G+0.1141B

式中g(i,j)表示的是转换后的灰度图像在点(i,j)处的灰度值，R、G、B为转换前的彩色图像在点(i,j)处的基色分量；

图像锐化，具体实现如下：梯度锐化：定义图像f(x,y)在点(x,y)处的梯度矢量G[f(x,y)]为一个二阶列矢量，梯度大的幅值即模值，即：

图像平滑处理，具体实现如下：

依次采用直方图均衡化的灰度增强算法和中值滤波法对图像进行平滑处理。

上述步骤4中的图像识别，灰度均值L的计算公式如下：

式中：M×N为图像的大小，f(i,j)为像素点(i,j)的灰度值，M为图像长度的像素点数目，N为图像宽度像素点数目，i＝0,1，2,3,4…M-1，j＝0,1，2,3,4…N-1。

上述步骤4中的图像识别，灰度共生矩阵为：

(1)角二阶矩：

式中：M₁为度量角二阶矩大小的参数，L为灰度均值，P(m,n,d,θ)为点对灰度值(i,j)出现的频率，a,b分别为灰度值j的像素点距离i的横坐标和纵坐标距离，θ是指灰度共生矩阵的生成方向，θ一般取0°、45°、90°、和135°，m＝1,2,3,4…L，n＝1，2,3,4…L；

(2)对比度：

式中：M₂为度量对比度大小的参数；

(3)相关性：

其中：

M₃为度量相关性大小的参数；

(4)熵：

式中：M₄为度量熵值大小的参数；

(5)逆差矩：逆差矩可以反映出图像纹理的局部变化，图像纹理越规则它的值越大；

式中：M₅为度量逆差矩大小的参数。

一种TBM破岩矿渣数字近景摄影测量识别及反演系统，包括TBM盾构掘进机破碎岩石产生的矿渣，排出矿渣所运送的输送带，三台分别设置在所述输送带中心位置水平线的左边、右边和正上方转台上的相机组成联机测量系统，以及根据映射关系建立的矿渣关系型数据库系统。不同矿渣图像的灰度均值和纹理特征不同，由BP神经网络模型学习并对矿渣进行识别分类。对分类的矿渣进行室内物理力学性质实验，得到矿渣对应的围岩岩性，单轴抗压强度σ_c、泊松比μ、变形模量E_b、裂隙发育程度分级C_r等数据。将场景重构后的矿渣图像的颜色信息(R，G，B)，纹理和室内实验得到的矿渣对应的围岩岩性，单轴抗压强度σ_c、泊松比μ、变形模量E_b、裂隙发育程度分级C_r建立映射关系，并由二者之间的联系建立矿渣关系型数据库，服务器主机设一旦遇到之前数据库没有遇到的矿渣，将会自动记录并人为进行数据更新。实时采集的目标检测图像，经过图像标定校正、图像采集、位置姿态解算、图像处理、图像识别分类、空间场景重构之后形成虚拟空间信息，在虚拟空间内对颜色、反射率和纹理等信息进行三维显示和采集测量后的要素以OSG图像对象形式显示在三维平台中，并对各项信息进行数值统计，生成并输出统计数据图表。将上述数据实时传至服务器主机，根据上述数据信息与矿渣数据库进行数据比对得到匹配度最高的各项特性参数。由TBM刀盘上掘进方向的总滚刀破岩力NF_n、护盾与地层间摩擦力F_f、后续设备与轨道间的摩擦力F_r从而建立硬岩TBM掘进推力预测模型；采用半理论半经验、多种计算方法结合的评价体系预测硬岩TBM掘进速度(PR)。将数据库匹配得到的各项特性参数包括R₁、σ_c、μ、E_b、C_r即时传给计算机，由TBM掘进推力和速度计算模型的结果对前方TBM掘进推力和速度做出预测，从而为推进系统设计及施工过程中参数实时调整提供重要依据。

本发明的优点在于采用矿渣数据库和现场数字近景摄影测量图像得到的矿渣各项参数进行比对，以此来推断岩石的物理力学参数，包括矿渣对应的围岩岩性，单轴抗压强度σ_c、泊松比μ、变形模量E_b、裂隙发育程度分级C_r等。继而利用TBM掘进推力和速度计算模型的计算结果对前方TBM掘进推力和速度做出预测，为推进系统设计及施工过程中参数实时调整提供重要依据，提高隧洞掘进开挖效率，降低施工成本。

附图说明

图1是本发明的流程框图。

图2是本发明TBM盾构掘进机施工示意图。

图3是输送带和三台相机位置的俯视图。

图4是两种不同岩石的灰度直方图。

图5是BP神经网络结构图。

图6是BP神经网络3层4结点网络结构图。

图中，切削器1；前盾构2；带状螺旋器器3；皮带4；联机摄像系统5；列车土箱6；矿渣7。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，下面结合附图对本发明的一个实施例作进一步的描述。

本发明所述的TBM破岩矿渣数字近景摄影测量识别及反演方法，包括以下步骤：

步骤1，图像标定校正

图像标定校正主要是由标定获得相机内参数矩阵和外参数矩阵，从而可以建立世界坐标系和摄像机坐标系以及图像坐标系的映射关系，再通过反变换以及灰度插值即可实现畸变图像的校正，标定步骤如下：

(1)在四个设站对室内高精度三维控制场进行摄影，用结构已知、高精度加工的标定尺作为空间参照物，使用直接线性算法(DLT)求解中心构想方程，求得相机内参数矩阵L。以求得的参数为初始值，考虑径向畸变和切向畸变，采用最小二乘法求解畸变参数k_c。图像坐标系和世界坐标系之间的转换公式如下：

令P＝K[R t]，得到：

式中：[R t]为相机外参数矩阵，包括旋转矩阵R和平移向量t，L称为相机内参数矩阵；

X_w、Y_w、Z_w为相机空间坐标；x，y为像点的像平面坐标；f为影像的内方位元素；

Z_c为像空间坐标，X,Y,Z为物方点的物方空间坐标；x₀，y₀为主点坐标；f_x，f_y分别为图像x轴和y轴上的尺度因子；γ是描述两图像坐标轴倾斜程度的参数；P＝K[R t]为3×4矩阵，通常称为相机矩阵或投影矩阵。

(2)建立相关的数学模型减少相机光学系统存在加工和装配误差，非线性畸变一般为几何畸变，即物点在像平面上的成像会偏离理想位置，典型的相机畸变模型可用下式表示：

式中：右边第一项是径向畸变，第二项是切向畸变；k₁,k₂,p₁,p₂分别是非线性畸变参数，与前面所述的成像模型中的L矩阵组成了非线性相机模型的内部参数。

(3)直接线性算法(DLT)是由共线方程式推导而来的直接建立坐标仪坐标与物方空间坐标间的关系式的一种算法，计算中不需要内外方位元素的初值。共线方程如下：

式中：x，y为像点的像平面坐标；(δ_x，δ_y)为物镜畸变；f为影像的内方位元素；X,Y,Z为物方点的物方空间坐标；X_s，Y_s，Z_s为摄站点的物方空间坐标；a_i，b_i，c_i(i＝1,2,3)为3个外方位角元素组成的9个方向余弦。由共线方程推导而来的DLT方程如下：

式中：l_i(i＝1，2，3……11)是相机内外方为元素的表达式；(δ_x，δ_y)是物镜畸变。该方程含11个未知数，故至少需6个控制点，才可解得l_i(i＝1，2，3……11)和畸变参数。然后再考虑旋转矩阵的正交性质，利用下式得到相机内参数矩阵L。

(4)试验时在四个设站对室内高精度三维控制场进行摄影，如图2所示。摄像装置的检校参数为：主距f像主点坐标(x₀,y₀)；畸变系数(k₁,k₂)为径向畸变系数；(P₁,P₂)为切向畸变系数。得到内参数矩阵L和畸变参数k_c；使用定标后的摄像装置内参数，求得的外方位元素，物方坐标真值依据共线条件方程推算像方点坐标；将实际量测的像方坐标经畸变改正后得到中误差像素大小。具体的畸变校正可以通过Matlab软件或者OpenCV实现。

步骤2，TBM盾构机破岩产生矿渣

如图2所示，盾构掘进机掘进时，首先使后盾构夹紧固定在隧道中，然后开动切削器1旋转，进行地层开挖，同时，开动前用盾构推进千斤顶，从而推动前盾构2前进。随着切削器1旋转开挖的同时，通过安装在切削器1周围的铲斗将矿渣7铲起，当铲斗内矿渣7通过带状螺旋器3旋转到顶部时，就自动掉入溜槽卸放到输送机上的皮带4，进而不断地排放到列车土箱6中送往地面。

步骤3，图像采集、姿态解算

步骤3.1，图像采集步骤为：如图3所示，将摄像装置固定在输送带中心位置水平线的左边、右边和正上方的转台上组成联机摄像系统，利用球铰可将相机调整一定范围角度，将照明装置与冲洗装置固定到合适位置。打开电源，通过摄像装置对输送带进行观察，如有必要可利用冲洗装置进行定向清洗，照明装置调整光照，至所成图像清晰可辨。矿渣在经过输送带4时，安装在输送带水平轴左边、右边和正上方的相机联机系统5通过控制测量和摄取输送带4上的矿渣影像，实时采集目标检测图像。

步骤3.2，姿态解算是指由一组2D点的映射中估计物体的3D姿态，其步骤如下：

相机的参数矩阵L已通过相机标定得到，本发明采用PNP算法进行空间姿态计算，其依据任意2D点的可视角度和对应3D点的角度一样。本文取4个控制点，即采用P4P算法。

一组2D和3D的对应点其中是单位方向矢量，通过标定矩阵L的逆将2D像素测量值变换为单位模长的3D方向矢量，得到：

未知量为从相机中心c到3D点p_i的距离d_i，其中由三角形Δ(c,p_i,p_j)余弦定理得:

式中：组合其它约束(f_ij，f_ik，f_jk)，并消去d_i,d_j可得如下关于的公式：

给定五个或者更多的对应点，可以产生(n-1)(n-2)/2个三元组来得到的值的线性估计(用SVD)。可以用连续两个估计的比率来估计，然后通过平均来获得d_i。确定d_i后，可得到一系列尺度自由的点然后通过3D点之间的刚体变换，可将其变换到对应的样本点{p_i}，从而完成线性投影矩阵确定。

步骤4，图像处理、识别及分类

步骤4.1，图像处理主要包括图像灰度化、图像锐化、图像平滑处理。其步骤包括：

1)为了提高矿渣图像的处理速度，减少计算量，需要对矿渣图像进行灰度化处理，图像灰度化就是指将彩色图像变成灰度图像的过程，RGB颜色模型中R(红)、G(绿)、B(蓝)三个基色分量决定了彩色图像像素的颜色，图像灰度化的方法如下所示：

R＝G＝B＝(W_RR+W_GG+W_BB)/3

其中W_B、W_R、W_G分别是R、G、B的权值。

2)因为人眼对红、绿、蓝这三个颜色的敏感度不同，当W_B＜W_R＜W_G时将会得到比较理想的灰度图像，常用的经验公式如下:

g(i,j)＝R＝G＝B＝0.2989+0.5870+0.1141

式中g(i,j)表示的是转换后的灰度图像在像素点(i,j)处的灰度值。

3)图像的锐化处理可以增强图像的边缘和图像的细节信息，梯度算子是基于一阶微分的图像增强。

梯度锐化：定义图像f(x,y)在点(x,y)处的梯度矢量G[f(x,y)]为一个二阶列矢量，梯度大的幅值即模值，即：

梯度的方向在f(x,y)最大变化率方向上，方向角可表示为：

在实际应用中，梯度的模可以用差分代替微分，如使用x,y方向上差分绝对值替代模来度量。

梯度的模(幅值)就是最大变化率方向的单位距离所增加的量。由梯度的计算可知，在图像灰度变化较大的边沿区域其梯度值大，在灰度变化平缓的区域梯度值较小，而在灰度均匀的区域其梯度值为零。我们根据得到的梯度值来返回像素的值，加强梯度值大的像素灰度值就可以突出细节了达到了锐化的目的。

4)如图4所示，本发明先采用直方图均衡化的灰度增强算法，后采用中值滤波法对图像进行平滑处理。

步骤4.2，图像识别的原理是不同的岩石在外观上有明显的不同，反映在不同的颜色、纹理等方面。图形特征提取主要包括岩石灰度均值和纹理的特征提取，从而实现矿渣岩石的图像识别。具体步骤如下:

1)设图像的大小为M×N，像素点(i,j)的灰度值为f(i,j)，灰度均值L的计算公式如下。

通过上式可以求出多幅煤图像的灰度均值，再利用这个公式求出多幅岩石图像的灰度均值，通过分析得到不同岩石图像灰度均值的范围。根据不同范围的图像灰度均值划分为不同种类的岩石。

2)本发明采取灰度共生矩阵进行纹理特征提取，灰度共生矩阵是像素距离和角度的矩阵函数，能反映出图像灰度在方向、相邻间隔、变化幅度上的综合信息。如果灰度共生矩阵的值集中在对角线附近，那么图像的纹理较为粗糙；如果灰度共生矩阵的值分散在矩阵各处，那么图像的纹理较细。

3)设任意一点像素(x,y)的灰度值为m和偏离它的一点(x+a,y+b)的灰度值为n两点组成一个点对，点对的灰度值即为(m,n)。使a和b保持不变，移动(x,y)就会得到不同的(m,n)。假设图像的灰度级为L,则(m,n)的组合就会有L×L种。计算出每一种(m,n)出现的频率P(m,n,d,θ)，灰度共生矩阵就是由每一种P(m,n,d,θ)构成的方阵，P(m,n,d,θ)的计算公式如下。

P(m,n,d,θ)＝{[(x,y),(x+a,y+b)|f(x,y)＝m；f(x+a,y+b)＝n}]

式中：θ是指灰度共生矩阵的生成方向，θ一般取0°、45°、90°、和135°。

4)将灰度共生矩阵进行归一化处理，公式如下：

P¹(m,n,d,θ)＝P(m,n,d,θ)/R

式中，P¹(m,n,d,θ)表示归一化后的灰度共生矩阵，R表示的是灰度共生矩阵中所有元素的累加和。然后再在这个基础上统计出特征参数。表示纹理的特征参数的主要有以下几种。

(1)角二阶矩：角二阶矩即能量，是用来对图像灰度分布的均匀性进行度量的特征参数，角二阶矩越大，从图像上看，纹理越粗糙；反之，纹理越细。

(2)对比度：对比度也叫主对角线惯性矩，它反应了图像局部灰度的变化情况，对比度越大，图像局部灰度变化越大，图像越清晰，纹理越明显，反之亦然。

(3)相关性：相关性是用来描述空间灰度共生矩阵的元素在行或列方向上的近似程度，纹理在某个方向越强，相关性越大，其值可以反映图像纹理的主方向。

(4)熵：熵可以反映图像纹理的复杂程度，当图像很光滑，没有纹理的时候，熵值为0；反之纹理越复杂，熵值越大。

(5)逆差矩：逆差矩可以反映出图像纹理的局部变化，图像纹理越规则它的值越大。

步骤4.3，图像分类：本发明采用BP神经网络结构的方法进行分类，因为一个隐层的BP网络可以逼近任何一个在闭区间内的连续函数，三层的BP神经网络就可以完成任意n维到m维的映射,所以选用三层神经网络结构，即输入层、中间层和输出层。BP神经网络结构图如图5所示。

为了方便说明，我们每层结点均取4，输入层的输入是处理过的GR、CNL、DEN和AC值，如图6所示。输出为花岗岩、泥岩、灰岩和砂岩的输出模式。网络分学习和测试两个阶段。具体过程如下：

1)网络初始化：给网络输入层、中间层和输出层各单元权向量U_ij、V_ij、W_ij赋初值，通常取[-1.0～1.0]之间的随机值。并为每层各单元阈值θ赋初值，θ可根据试验定，通常取[1.0～5.0]。设置学习速度为μ，取[0～1.0]。

2)网络提供的样本模式对X＝(X₀～X₃)，对应的期望值D＝(d₀,d₁,d₂,d₃)和要求学习的次数。

3)输入层各单元的输出值Y_i为：

式中：i＝1～3，f为S型函数f(x)＝(1+e^-x)^-1

中间各结点的输出为：

输出层各结点的输出为：

4)计算输出层、中间层和输入层的权向量，如果d_i和Out_i分别表示输出层的第m个结点的期望值和实际值，Z_i，Y_i分别表示中间层和输入层的结点的实际输出，e₃、e₂、e₁分别表示输出层、中间层和输入层的误差，则权项调节值为：

输出层：ΔW_ij＝μe_3iZ_j i＝0～3，j＝0～3

中间层：ΔV_ij＝μe_2iZ_j i＝0～3，j＝0～3

输入层：ΔU_ij＝μe_1iZ_j i＝0～3，j＝0～3

其中μ为学习效率。

5)将所提供的样本作输入进行样本训练学习，直到满足收敛条件为止。然后从第(3)步开始，继续输入新的样本，如簇循环，直到输入足够的样本为止，并存储联接矩阵系数。

6)利用所存储的联接矩阵对新的样本进行测试识别。

不同矿渣的图像灰度均值L大小和纹理特征不同，由不断训练的BP神经网络结构对矿渣进行分类。

步骤5，对已分类的矿渣进行室内材料与力学性质实验，得到矿渣对应的围岩岩性类别R₁、单轴抗压强度σ_c、泊松比μ、变形模量E_b、裂隙发育程度分级C_r。

步骤6，将场景重构后的矿渣图像的灰度均值，纹理特征和室内实验得到的矿渣对应的围岩岩性类别R₁、单轴抗压强度σ_c、泊松比μ、变形模量E_b、裂隙发育程度分级C_r建立映射关系，并由二者之间的联系建立的矿渣关系型数据库。

步骤7，在虚拟空间内对颜色、反射率和纹理等信息进行三维显示和采集测量后的要素以OSG图像对象形式显示在三维平台中，并对各项信息进行数值统计，生成并输出统计数据图表。三维平台在输出得到矿渣表面的颜色和纹理等信息图表后，将上述数据实时传至服务器主机，所述主机在矿渣数据库中搜索与之信息最接近的岩石或土样，从而根据信息对应关系推测前方岩石的参数，包括围岩岩性类别R₁、单轴抗压强度σ_c、泊松比μ、变形模量E_b、裂隙发育程度分级C_r等。服务器主机一旦遇到之前数据库没有遇到的矿渣，即主机在矿渣数据库中搜索不到与之信息接近的岩石或土样，会将新的图像信息即矿渣表面的颜色和纹理等信息图表数据加入矿渣关系型数据库中，更新矿渣关系型数据库。

步骤8，在硬岩TBM掘进过程中，快速准确地计算硬岩TBM在掘进过程中所需要的总推力值和预测TBM的掘进速度，可以为推进系统设计及施工过程中参数实时调整提供重要依据。

本发明的硬岩TBM掘进过程中预测总推力的计算方法包括以下步骤：

(1)计算TBM刀盘上掘进方向的滚刀破岩力F_n：

式中，r为滚刀半径；δ为滚刀与岩石的接触角，有δ＝arccos(r-d)/r)，ω为滚刀刀圈宽度，p₀为位于刀具正下方的基准压力；ψ为接触压力分布指数，取决于滚刀和岩石的参数。

通过一系列实验测得的各种工况不同岩石滚刀正下方的的基准压力，由统计回归分析后得到了基准压力的计算公式：

式中，δ是滚刀岩石接触角，S是刀间距，σ_c是岩石单抗压强度，σ_i是岩石抗剪强度，C≈2.12，是无量纲系数，ω为刀圈宽度。

接触压力的分布形态会随滚刀的贯入深度发生变化，因而可以根据滚刀贯入深度的变化调整ψ的取值，提高计算结果的准确性。本发明采用用最小二乘拟合思想，认为只要使计算值与试验值相对误差的平方和取极小值，即获得了最优ψ值，然后对此时的ψ进行回归分析，便得到了以贯入深度d为自变量的ψ的回归计算公式。

TBM刀盘上通常装有数十把滚刀，各滚刀交错排列。尽管同一时刻不同位置滚刀受力略有差别，但从统计角度分析刀盘整体受力时，可认为各把滚刀平均受力近似相等，那么刀盘整体破岩力可近似表示为：

式中，N为刀盘上装有的滚刀数量。

(2)计算护盾与地层间摩擦力F_f：

式中，μ₁是护盾与岩石的摩擦系数；D是刀盘直径(m)；l是护盾长度(m)；是护盾与洞壁接触范围(rad)；p是接触压力(Mpa)，由护盾与地层接触情况决定。

(3)计算后续设备与轨道间的摩擦力F_r：F_r＝μ₂mg

式中，μ₂是后续设备与轨道间的摩擦系数；m是TBM后续设备质量(单位：t)，g是重力加速度，一般取9.8m/s²。

(4)计算TBM掘进总推力F:

式中，N为刀盘上装有的滚刀数量，π为圆周率，ω为滚刀刀圈宽度，r为滚刀半径，d为贯入深度(m)，p₀为位于刀具正下方的基准压力；n＝0,1,2,…,∞；ψ为接触压力分布指数，取决于滚刀和岩石的参数；δ是滚刀岩石接触角；D是刀盘直径(m)；l是护盾长度(m)；是护盾与洞壁接触范围(rad)；p是接触压力(Mpa)，由护盾与地层接触情况决定；μ₁是护盾与岩石的摩擦系数，μ₂是后续设备与轨道间的摩擦系数；m是TBM后续设备质量(单位：t)，g是重力加速度，一般取9.8m/s²。

(5)计算TBM掘进速度：

根据数量化理论M的原理进行数据预处理，采用Spss多元线性回归求解工具，以最小二乘法理论为依据，在初始设置解释变量的基础上，逐步删除不显著的解释变量，求解正则方程，最终得到定性参数包括围岩岩性类别R₁等1个解释变量，定量参数包括单轴抗压强度σ_c、泊松比μ、变形模量E_b、裂隙发育程度分级C_r等4个解释变量的预测模型，正则方程如下：

Y(PR)＝β₀+β₁(σ_c)+β₂(R₁)+β₃(μ)+β₄(E_b)+β₅(C_r)

其中，Y(PR)为因变量，在这里为掘进速度；β₀、β₁、β₂、β₃、β₄、β₅为各项自变量对应的系数。根据收集的现场试验数据，代入上述的正则方程，经过多元线性回归计算可以得到模型参数估计结果，得到的系数矩阵为：

求解得到自变量系数越大，则说明系数的影响因子越大，反之则越小；自变量系数是正的，则说明该自变量对TBM掘进速度呈显著正效应，若是负的，则呈负效应。

由多元线性回归分析得到的各项自变量系数后，经过图像处理和识别分类的矿渣三维图像信息关系型数据库进行数据比对得到匹配度最高的各项特性参数，将各项特性参数包括R₁、σ_c、μ、E_b、C_r代入TBM掘进速度预测模型中即可得到掘进速度Y(PR)。

将最新得到的R₁、σ_c、μ、E_b、C_r等现场试验数据加入之前的数据表，从而不断通过多元线性回归修正并更新TBM掘进速度预测模型，并且将所述模型用于下一回合TBM掘进速度的预测，不断进行数据回归、预测速度、更新模型的循环。

现场的施工技术员可依据以上三种TBM掘进速度计算模型的结果对前方TBM掘进速度做出预测。

应当理解的是，对本领域普通技术人员来说，可以根据上述说明加以改进或变换，而所有这些改进和变换都应属于本发明所附权利要求的保护范围。

Claims (10)

1.一种TBM破岩矿渣数字近景摄影测量识别及反演方法，其特征在于，包括以下步骤：

步骤1，图像标定校正：主要是由标定获得相机内参数矩阵和外参数矩阵，从而建立世界坐标系和摄像机坐标系以及图像坐标系的映射关系，再通过反变换以及灰度插值实现畸变图像的校正；

步骤2，TBM盾构机破岩产生矿渣；

步骤3，图像采集、姿态解算；

图像采集：在TBM破岩矿渣输送带的左边、右边和正上方的转台上分别固定摄像装置，所述摄像装置组成联机摄像系统，从输送带的左、右、上方采集矿渣图像；

姿态解算：采用PNP算法进行空间姿态计算，依据任意2D点的可视角度和对应3D点的角度一样，由一组2D点的映射中估计物体的3D姿态；

步骤4，图像处理、识别及分类；

图像处理：包括对矿渣图像进行灰度化处理，将彩色矿渣图像变成灰度矿渣图像；

图像识别：对灰度矿渣图像进行岩石图像灰度均值和纹理特征提取；计算多幅岩石图像的灰度均值，通过分析得到不同岩石图像灰度均值的范围，根据不同范围的图像灰度均值划分为不同种类的岩石；采用灰度共生矩阵进行纹理特征提取；

图像分类：采用不断训练的BP神经网络结构对矿渣进行分类；岩石灰度均值、纹理特征是作为BP神经网络结构的输入值数据，用来训练区分不同的岩石图像；图像分类的处理对象是岩石图像的灰度均值和纹理特征等提取的图像信息，处理结果是根据不同的图像信息分类矿渣图像，从而识别对应的矿渣；

步骤5，通过对矿渣图像的分类区分出所对应的矿渣种类，对已分类的矿渣进行室内物理力学实验，得到矿渣对应的围岩岩性特征参数，包括单轴抗压强度σ_c、围岩岩性类别R₁、泊松比μ、变形模量E_b、裂隙发育程度分级C_r；

步骤6，建立图像信息和特征参数之间的矿渣关系型数据库；所述图像信息是指矿渣图像的灰度均值和纹理特征的数字化信息，所述特征参数是指矿渣对应的围岩岩性特征参数；

步骤7，在矿渣关系型数据库搜索相似的图像信息，若存在相似图像信息，将矿渣参数上传至服务器主机储存的关系型数据库中，并据此推测现场矿渣对应的围岩岩性特征参数；若不存在相似图像信息，将新的图像信息加入矿渣关系型数据库中，更新矿渣关系型数据库；

步骤8，建立TBM掘进推力和掘进速度预测模型，预测前方TBM掘进推力和速度；

TBM掘进总推力F预测模型如下：

式中，N为刀盘上装有的滚刀数量；π为圆周率；ω为滚刀刀圈宽度；r为滚刀半径；d为贯入深度；p₀为位于刀具正下方的基准压力；n＝0,1,2,…,∞；ψ为接触压力分布指数，取决于滚刀和岩石的参数；δ是滚刀岩石接触角；D是刀盘直径；l是护盾长度；是护盾与洞壁接触范围；p是接触压力，由护盾与地层接触情况决定；μ₁是护盾与岩石的摩擦系数，μ₂是后续设备与轨道间的摩擦系数；m是TBM后续设备质量，g是重力加速度；

TBM掘进速度预测模型如下：

Y(PR)＝β₀+β₁(σ_c)+β₂(R₁)+β₃(μ)+β₄(E_b)+β₅(C_r)

其中，Y(PR)为因变量，为掘进速度；β₀、β₁、β₂、β₃、β₄、β₅为各项自变量对应的系数；根据收集的现场试验数据，代入上述的正则方程，经过多元线性回归计算得到模型参数估计结果，得到的系数矩阵为：

将步骤7中推测得到的现场矿渣对应的围岩岩性特征参数单轴抗压强度σ_c、围岩岩性类别R₁、泊松比μ、变形模量E_b、裂隙发育程度分级C_r代入上述TBM掘进速度计算模型中即得到掘进速度Y(PR)。

2.根据权利要求1所述的TBM破岩矿渣数字近景摄影测量识别及反演方法，其特征在于：

上述步骤8中，将最新得到的R₁、σ_c、μ、E_b、C_r现场试验数据加入之前的数据表，从而不断通过多元线性回归修正并更新TBM掘进速度计算模型，并且将所述模型用于下一回合TBM掘进速度的预测，不断进行数据回归、预测速度、更新模型的循环。

3.根据权利要求1所述的TBM破岩矿渣数字近景摄影测量识别及反演方法，其特征在于：

上述步骤1中，图像坐标系和世界坐标系之间的转换公式如下：

其中，P＝K[R t]为相机矩阵或投影矩阵；[R t]为相机外参数矩阵，包括旋转矩阵R和平移向量t；L为相机内参数矩阵，x₀，y₀为主点坐标；f_x，f_y分别为图像x轴和y轴上的尺度因子，γ是描述两图像坐标轴倾斜程度的参数；Z_c为像空间坐标；x，y为像点的像平面坐标；X,Y,Z为物方点的物方空间坐标。

4.根据权利要求3所述的TBM破岩矿渣数字近景摄影测量识别及反演方法，其特征在于：

上述步骤1中，x₀，y₀，f_x，f_y计算公式如下：

式中：x₀，y₀为主点坐标；f_x，f_y分别为图像x轴和y轴上的尺度因子；l_i是相机内外方为元素的表达式，i＝1，2，3……11。

5.根据权利要求3所述的TBM破岩矿渣数字近景摄影测量识别及反演方法，其特征在于：

上述步骤1中，摄像装置的检校参数为：主距f像主点坐标(x₀,y₀)；畸变系数(k₁,k₂)为径向畸变系数；(P₁,P₂)为切向畸变系数；得到内参数矩阵L和畸变参数k_c；使用定标后的摄像装置内参数，求得的外方位元素，物方坐标真值依据共线条件方程推算像方点坐标；将实际量测的像方坐标经畸变改正后得到中误差像素大小。

6.根据权利要求1所述的TBM破岩矿渣数字近景摄影测量识别及反演方法，其特征在于：

上述步骤3的图像采集步骤，设置球铰调整摄像装置摄像角度，设置照明装置提供光照使所成图像清晰可辨，设置冲洗装置对摄像头进行清洗。

7.根据权利要求1所述的TBM破岩矿渣数字近景摄影测量识别及反演方法，其特征在于：

上述步骤3的姿态解算步骤，具体实现如下：

通过相机标定来获取相机内参数即基础矩阵和畸变向量，已知物体三维空间坐标系的点以及对应的图像像素坐标点；相机的参数矩阵L已通过相机标定得到，故可用PNP算法进行空间姿态计算，其依据任意2D点的可视角度和对应3D点的角度一样；采用P4P算法，取4个控制点：

一组2D和3D的对应点其中是单位方向矢量，通过标定矩阵L的逆将2D像素测量值变换为单位模长的3D方向矢量，得到：

未知量为从相机中心c到3D点p_i的距离d_i，其中由三角形Δ(c,p_i,p_j)余弦定理得：

式中：组合其它约束(f_ij，f_ik，f_jk)，并消去d_i,d_j可得如下关于的公式：

给定n个对应点，可以产生(n-1)(n-2)/2个三元组来得到的值的线性估计；可以用连续两个估计的比率来估计，然后通过平均来获得d_i；确定d_i后，可得到一系列尺度自由的点然后通过3D点之间的刚体变换，可将其变换到对应的样本点{p_i}，从而完成线性投影矩阵确定。

8.根据权利要求1所述的TBM破岩矿渣数字近景摄影测量识别及反演方法，其特征在于：

上述步骤4中的图像处理，在图像灰度化处理后，还进行了图像锐化和图像平滑处理：

图像灰度化，具体实现如下：

g(i,j)＝R＝G＝B＝0.2989R+0.5870G+0.1141B

式中g(i,j)表示的是转换后的灰度图像在点(i,j)处的灰度值，R、G、B为转换前的彩色图像在点(i,j)处的基色分量；

图像锐化，具体实现如下：梯度锐化：定义图像f(x,y)在点(x,y)处的梯度矢量G[f(x,y)]为一个二阶列矢量，梯度大的幅值即模值，即：

图像平滑处理，具体实现如下：

依次采用直方图均衡化的灰度增强算法和中值滤波法对图像进行平滑处理。

9.根据权利要求1所述的TBM破岩矿渣数字近景摄影测量识别及反演方法，其特征在于：

上述步骤4中的图像识别，灰度均值L的计算公式如下：

式中：M×N为图像的大小，f(i,j)为像素点(i,j)的灰度值，M为图像长度的像素点数目，N为图像宽度像素点数目，i＝0,1，2,3,4…M-1，j＝0,1，2,3,4…N-1。

10.根据权利要求1所述的TBM破岩矿渣数字近景摄影测量识别及反演方法，其特征在于：

上述步骤4中的图像识别，灰度共生矩阵为：

(1)角二阶矩：

式中：M₁为度量角二阶矩大小的参数，L为灰度均值，P(m,n,d,θ)为点对灰度值(i,j)出现的频率，a,b分别为灰度值j的像素点距离i的横坐标和纵坐标距离，θ是指灰度共生矩阵的生成方向，θ一般取0°、45°、90°、和135°，m＝1,2,3,4…L，n＝1，2,3,4…L；

(2)对比度：

式中：M₂为度量对比度大小的参数；

(3)相关性：

其中：

M₃为度量相关性大小的参数；

(4)熵：

式中：M₄为度量熵值大小的参数；

(5)逆差矩：逆差矩可以反映出图像纹理的局部变化，图像纹理越规则它的值越大；

式中：M₅为度量逆差矩大小的参数。