本发明涉及高速铁路隧道围岩等级划分方法，具体涉及一种多井多参数约束下隧道波速反演成像及围岩分级的方法。
背景技术：
近年来，随着高速铁路的发展，山区铁路的建设日益增多，相应的铁路长大埋深隧道日益增多，而该类隧道的围岩完整性是当今世界各国工程界亟待解决的问题。
在隧道设计阶段，围岩完整性评价主要通过地质调查法、岩芯钻探测试法、岩体纵波波速法相结合起来综合确定。
地质调查是通过收集区域地质资料并进行现场踏勘，分析线路通过位置的地质构造、地质年代、地层岩性进行隧道围岩等级划分。该方法仅能看到地表的岩性特征，无法获得深部围岩的情况，只能根据经验推断深部洞身的围岩等级，而且还无法考虑地下水的情况。
利用钻孔资料进行围岩等级划分效果较好，可以根据钻孔取得岩芯直观的分析岩石的完整性及强度，从而给定围岩等级。但是对于长大埋深隧道来说却是困难的。首先，长大埋深隧道一般地形条件比较复杂，线路通过位置很多都是人迹罕至地区，搬运钻机上山需要修路，修路的周期很长；其次，既然是深埋，钻孔成孔周期很长，深孔往往后滞整个地质勘察，最后，即使取得了钻孔资料，但毕竟仅仅式一个小孔的岩石特征，更广大的、区域的岩性特征不能体现出来。
在铁路隧道设计阶段，长大埋深隧道的围岩纵波波速资料一般很难获取，目前获取弹性波波速的主要方法有：地震勘探、测井。
隧道一般都是跨越山脉，受地形条件的制约，地震方法很难开展，尤其是在深长隧道全面开展地震勘探工作，不仅成本高，而且难度很大。另外，在复杂地区，尤其是高陡地区，地震资料必然会受到影响，所以获得的速度数据精度受到影响。同时由于国家对易燃易爆物品的管制和爆炸震源的安全隐患，也限制了深部弹性波资料的获取。
测井方法可以准确的获取地层的速度数据。但是测井获得的速度数据往往只是一个点位上的数据，无法获取整条隧道的围岩速度场。如果加密钻孔，不仅成本高，而且会加大野外工作量，拖长工期。
如何利用电性参数进行围岩等级划分化，受到越来越多关注。希瑞研究利用电性参数与岩石完整性关系，利用隧道剖面电阻率参数判断隧道围岩岩石完整性进行围岩等级划分；许广春研究了利用电性差异进行围岩等级划分。但目前为止，以电性参数进行围岩等级划分尚无统一的量化标准。
综上所述，在勘察设计阶段，弹性波波速的获取亟须解决。
技术实现要素：
本发明的目的是提供一种多井多参数约束下隧道波速反演成像及围岩分级的方法，解决勘察设计阶段铁路隧道中难以获取隧道围岩波速信息的问题，为设计阶段围岩等级划分提供数据支撑，为前期方案研究、方案比选论证提供依据。
本发明所采用的技术方案为：
多井多参数约束下隧道波速反演成像及围岩分级的方法，其特征在于：
包括以下步骤：
第一步：多井多参数约束下的隧道波速反演成像：
利用隧道电性参数资料在控制钻孔岩芯标本物性参数的约束下，通过反演得到隧道围岩波速参数，实现隧道围岩波速成像；
第二步：多井多参数约束下的围岩分级：
利用控制钻孔岩芯完整性与波速参数的对应关系，建立隧道围岩波速分级标准，在隧道波速成像的基础上，实现隧道围岩的围岩等级划分。
第一步，多井多参数约束下的隧道波速反演成像，由以下步骤实现：
步骤1：利用faust公式
v≈kρahβ(1)
式中k为常数；α,β为权重；ρ为岩石电阻率；h为岩石埋深；v为岩石波速；将电阻率参数求取对数的方式进行量纲化处理，即用lgρ代替(1)式中的ρ参数；
将(1)式用多项式展开，并舍弃高次项，得(2)式：
v≈a1(lgρ)hβ+a2((lgρ)hβ)2+a3((lgρ)hβ)3+a4((lgρ)hβ)4+a5(2)
步骤2：将m个控制钻孔电阻率、深度及对应的波速值分别代入(2)式，得到一组关于系数(aj1,aj2,aj3,aj4,aj5)j＝1,m的超定方程组：
式中
式中m为钻孔个数，n为各井中选取约束条件个数；
写成矩阵形式为v＝fa(4)
a＝(aj1aj2aj3lajn)j＝1,2,l,m
其中v为波速值矩阵；f为电性参数矩阵；a为待求系数矩阵；
采用自适应加权最小二乘法能对数据进行滤波，自适应加权算法是在矩阵中引入自适应加权函数；从而减小扰动数据对方程组的影响；
选择权函数：
式中ρci,hci为钻孔控制点位电阻率与深度；ρi,hi为计算点位电阻率与深度；(ρcihci)rms为钻孔控制点位电阻率与深度乘积的均方根值；
在(4)式中左乘权系数矩阵，引入自适应加权函数，依据矩阵运算法则，改为
wv＝wfa(5)
步骤4：利用奇异值分解法求解(5)式，得各控制钻孔中波速反演系数；
将深度加权系数、反演系数和控制钻孔中岩芯电性参数代入式(2)，获得井中电性参数反演波速成果；
步骤5：将反演得到的波速vj与井中实测的波速vp比较，求出二者的方差：
式中vp为井中实测波速参数、vj为反演波速参数，n为井中约束条件个数；
将深度加权系数矩阵[β1lβm]按台劳级数展开，并忽略二次项以上的项，式(6)表达改为求深度加权系数修改量；使目标函数趋于极小：
依据多元函数极值原理，上式要趋于极小,则对于各参数要满足下面的线性方程:
用差分方法来计算偏导数
取,δβj＝0.1βj则βj的偏导数为
vj(β1,,..βj....,,βm)为原反演系数和深度加权系数反演波速，
精度不满足要求时，求取修改量后利用式(10)获得新的深度加权系数返回步骤2，重新进行反演；
另外为了防止参数修改过量，实际过程中作如下规定:δβj＞1.2βj时取δβj＝1.2βj；δβj＜-0.2βj时取δβj＝-0.2βj，即每次修改量不超过原有模型参数值的0.2倍,保证收敛稳定；
步骤6：利用步骤5反复迭代深度加权系数，使反演波速逐步逼近井中实测波速；当精度满足要求，求得各控制钻孔的反演系数和深度加权系数；
分段计算算法将按井求取系数矩阵，获得局部的最优化参数，然后采用“就近原则”即，该电性参数里程坐标距某个控制钻孔距离较近，就采用该钻孔的反演系数及深度加权系数；
假设控制钻孔zk1,zk2,l,zkm，里程坐标分别为x1,x2,l,xm，x为电性参数里程坐标；若|x-xi|＝min，采用zki控制钻孔的反演系数及深度加权系数；
将深度加权系数、反演系数和隧道剖面电性参数代入式(2)，获得隧道剖面电性参数反演波速剖面成果；
步骤7：通过visualc++对surfer、autocad进行二次开发，实现隧道围岩波速反演成像。
第二步：多井多参数约束下的围岩分级，，由以下步骤实现：
步骤1：采集控制钻孔中各岩性的岩芯标本，并记录岩芯标本的取样深度；
步骤2：在室内测试岩芯标本的波速参数及电阻率参数，并统计各岩性波速参数及电阻率参数；依据标本的完整性及对应的波速参数，建立隧道围岩的波速分级标准；
步骤3：依据铁路工程物探规范、井中岩石完整性与波速对应关系及地质资料，完成隧道围岩等级的划分。
本发明具有以下优点：
本发明涉及的技术方法，是利用电性参数、井中参数反演剖面波速参数进行围岩等级划分的一种方法，能够利用少量的控制钻孔及剖面电性资料获得隧道剖面上波速信息，有效的解决勘察设计阶段铁路隧道中难以获取隧道围岩波速信息的问题。
附图说明
图1波速反演成像成果图。
图2围岩等级划分成果图。蓝色实线为控制钻孔位置及深度，红色实线为隧道洞身位置，罗马数字为围岩等级。
图3分段计算算法示意图。
具体实施方式
下面结合具体实施方式对本发明进行详细的说明。
本发明涉及一种多井多参数约束下隧道波速反演成像及围岩分级的方法，包括多井多参数约束下的隧道波速反演成像及多井多参数约束下的围岩分级方法两个部分。
第一部分利用数学推导，建立控制钻孔建立隧道围岩电性参数与波速参数的数学关系，在钻孔岩芯的约束下，实现波速反演成像(如图1所示)。反演结果获得了隧道剖面上连续的波速信息，解决了隧道围岩波速难以获取的问题，为设计阶段围岩等级划分提供数据支撑。
第二部分采集控制钻孔中岩芯标本，统计隧道围岩各岩性的波速参数，并通过标本的完整程度及相应波速参数，建立隧道围岩的波速分级标准，最终利用波速反演成像结果，实现隧道围岩的等级划分(如图2所示)。围岩划分结果为前期方案研究、方案比选论证提供依据。
第一部分：多井多参数约束下的隧道波速反演成像方法
下面结合具体实施方式对本发明进行详细的说明。
步骤1：利用faust公式v≈kρahβ(1)
式中k为常数；α,β为权重；ρ为岩石电阻率；h为岩石埋深；v为岩石波速。将电阻率参数求取对数进行量纲化处理，即用lgρ代替(1)式中的ρ参数。
将(1)式用多项式展开，并忽略高次项，可得到(2)式：
v≈a1(lgρ)hβ+a2((lgρ)hβ)2+a3((lgρ)hβ)3+a4((lgρ)hβ)4+a5(2)
步骤2：将m个控制钻孔电阻率、深度及对应的波速分别代入(2)式，即可得到一组关于系数(aj1,aj2,aj3,aj4,aj5)j＝1,m的超定方程组：
式中
其中m为钻孔个数，n为各钻孔选取约束条件个数。
写成矩阵形式为v＝fa(4)
其中v为波速值矩阵；f为电性参数矩阵；a为待求系数矩阵。
步骤3：采用自适应加权最小二乘法能对数据进行滤波，自适应加权算法是在矩阵中引入自适应加权函数。从而减小扰动数据对方程组的影响。
自适应加权函数的选择原则：按照其偏离均方根值的距离，对数据进行加权，偏离距离值越大，其权系数越小，对方程组的“贡献”越小；偏离距离值越小，其权系数越大，对方程组的“贡献”越大；为非负的关于原点对称的函数。
此次选择的权函数为：
式中ρci,hci为钻孔控制点位电阻率与深度；ρi,hi为控制点位电阻率与深度；(ρcihci)rms为钻孔控制点位电阻率与深度均方根值。
在(4)式中左乘权系数矩阵，引入自适应加权函数，依据矩阵运算法则，可改为
wv＝wfa(5)
步骤4：利用奇异值分解法求解(5)式，可得各控制钻孔中波速反演系数。
奇异值分解法认为任意n×m阶矩阵f均可分解为f＝uwvt,这里u为n×m阶正交阵和v为m×m阶正交阵。
其中δ1,l,δr为矩阵f的奇异值,r是矩阵的秩.当f非奇异时,奇异值较大,方程组(5)有广义逆解,f＝vw-1ut这里
当f接近奇异时,有的奇异值就较小,此w-1时由于系数过大,上述解的误差就较大。为了解决这个问题,维根斯(wiggins)提出用最接近的矩阵r来代替f,,而r＝uwevt
w中小的奇异值在这里便被零代替了，因此有较精确的广义逆解：
将深度加权系数、反演系数和控制钻孔中岩芯电性参数代入式(2)，即可获得控制钻孔中电性参数反演波速成果。
步骤5：将反演得到的波速vj与井中实测的波速vp比较，求出二者的方差：
式中vp为井中实测波速参数、vj为反演波速参数，n为井中约束条件个数。
将深度加权系数矩阵[β1lβm]按台劳级数展开，并忽略二次项以上的项，式(6)表达可改为求深度加权系数修改量。使目标函数趋于极小：
依据多元函数极值原理，上式要趋于极小,则对于各参数要满足下面的线性方程:
可用差分方法来计算偏导数
取,δβj＝0.1βj则βj的偏导数为
vj(β1,,..βj....,,βm)为原反演系数和深度加权系数反演波速，
利用松弛迭代算法方法求解得到深度加权系数修改量。
将δa{δaj,j＝1,2,l,5}分别代入式(14)可得反演参数修改量的方程组
根据sor原理，将(14)式改写为：
(15)式为sor迭代公式，式中w为松弛因子；迭代的收敛条件为
当w<1，为低松弛迭代；
当w＝1，为松弛迭代；
当w>1，为超松弛迭代。
精度不满足要求时，求取修改量后利用式(15)获得新的深度加权系数返回步骤2，重新进行反演。
另外为了防止参数修改过量，实际过程中作如下规定:δβj＞1.2βj时取δβj＝1.2βj；δβj＜-0.2βj时取δβj＝-0.2βj，即每次修改量不超过原有模型参数值的0.2倍,保证收敛稳定。
步骤6：利用步骤5反复迭代深度加权系数，使反演波速逐步逼近控制钻孔岩芯标本波速；当精度满足要求，即可求得各控制钻孔的反演系数和深度加权系数。
分段计算算法将按井求取系数矩阵，获得局部的最优化参数，然后采用“就近原则”，该电性参数里程坐标距某个控制钻孔井距离较近，就采用该井的反演系数及深度加权系数。
图3中zk1,zk2,l,zkm为控制钻孔，里程坐标分别为x1,x2,l,xm，小圆点为某个电性剖面测点，x为电性参数里程坐标。若|x-xi|＝min，采用zki控制钻孔的反演系数及深度加权系数。
将深度加权系数、反演系数和隧道剖面电性参数代入式(2)，即可获得隧道剖面电性参数反演波速剖面成果。
步骤7：通过visualc++对surfer、autocad进行二次开发，实现隧道围岩波速反演成像。
第二部分：多井多参数约束下的围岩分级方法
步骤1：采集控制钻孔中各岩性的岩芯标本，并记录岩芯标本的取样深度。
步骤2：在室内测试岩芯标本的波速参数及电阻率参数，并统计各岩性波速参数及电阻率参数。依据标本的完整性及对应的波速参数，建立隧道围岩的波速分级标准。
步骤3：依据铁路工程物探规范、井中岩石完整性与波速对应关系及地质资料，完成隧道围岩等级的划分。
本发明的内容不限于实施例所列举，本领域普通技术人员通过阅读本发明说明书而对本发明技术方案采取的任何等效的变换，均为本发明的权利要求所涵盖。