CN104407394A - 基于电阻率跨孔ct的地铁盾构区间孤石精细化探测方法 - Google Patents

Abstract

本发明公开了一种基于电阻率跨孔CT的地铁盾构区间孤石精细化探测方法，包括：在掘进面前方未开挖段设计科学合理的钻孔平面布置方案；铺设电阻率跨孔CT测线网，采用二维探测普查、三维探测详查相结合的探测方案；数据采集利用新型组合式观测模式，获取更多关于孤石异常体的有效信息；通过构造携带距离加权函数的三维电阻率反演目标函数及反演方程，采用偏导数矩阵的并行解析快速求解算法；形成测区范围内孤石异常体探测的电阻率成像结果图，并结合已有地质分析推断孤石体的大小、空间分布等情况。本发明可实现孤石体异常体的三维高分辨率精确成像，为地铁盾构掘进机施工隧道安全快速掘进提供了有力的技术支持。

Description

基于电阻率跨孔CT的地铁盾构区间孤石精细化探测方法

技术领域

本发明涉及一种地铁盾构区间孤石精细化探测方法，尤其涉及一种基于电阻率跨孔CT的、在地铁盾构掘进机施工隧道条件下对隧道掘进面前方孤石体或孤石群进行精细化探测的方法。

背景技术

近年来，我国城市地铁工程建设迎来了宝贵的机遇期和建设高潮。在地铁施工过程中，盾构掘进机施工方法逐渐成为我国地铁工程建设的重要趋势。盾构掘进机施工具有“掘进速度快、施工扰动小、成洞质量高、综合经济社会效益高、施工安全文明”的优点，但其对地层的适应性较差。地铁工程勘察中发现，盾构掘进机在长距离掘进过程中经常遇到“孤石群地层”，特别是在我国广州、深圳、长沙地铁施工中较为常见。未探明的孤石群会给地铁盾构施工带来重大安全隐患，不仅导致盾构刀盘频繁被卡甚至严重变形，延误工期造成巨大经济损失，严重时还会导致工作面喷涌、塌方，引起突发地质灾害。因此，在盾构掘进前，尽可能准确地探查清楚掘进区域地铁隧洞洞身范围内孤石群的赋存状态，是地铁盾构工程施工过程中迫切需要解决的难题。

钻探和物探方法是探测孤石的重要手段。钻探可以直观地揭露地铁盾构区间局部点位的地下球状风化体的竖向大小、赋存位置、表现形式等信息，但对于连续区域的探查则具有“一孔之见”的局限性，通过钻探能揭露到的孤石是十分有限的，即使通过加密钻孔提高揭露孤石的机率，但受成本、场地条件等限制，是不现实的。

目前，用于孤石探查的物探方法主要有：地质雷达法、地震反射波法、跨孔电阻率CT法等，由于盾构施工方法和装备所限，在洞内实施超前探测很难具备可用的工作空间和有效的观测条件，因此目前主要以地表探测手段为主。采用地面物探方法探测孤石，不同方法的探测结果存在较大差异，均达不到理想效果。对于有一定深度的目标体，孔间/孔中物探方法具有相对较好的探测效果，其分辨率远高于地面装置型式。电阻率跨孔CT方法可以采集丰富的数据量，具有高分辨率精细探查的显著优势，在孤石定位方面取得了一定的探测效果。

对于利用电阻率CT法探测地铁孤石而言，分辨率、定位精度和可靠性是衡量其探测效果的关键指标。由于在电阻率CT探测孤石的多解性压制方法、有效观测方式及探测效果影响因素等方面缺乏系统深入的研究，导致缺乏一套电阻率CT探测孤石的合理科学的探测方案和工作标准，时常出现“探不到”(即难以识别小尺寸孤石)、“探不准”(即对孤石难以准确定位)等问题，使得孤石探测的可靠性难以保证，无法满足工程实际需要。具体而言，应用电阻率CT探测孤石存在以下几个关键性问题：

①同其他地球物理方法一样，电阻率CT探测方法的反演与解译存在着多解性难题，严重时会导致反演结果与实际情况差别较大，对孤石难以准确定位，因此亟需改善跨孔电阻率CT反演问题的多解性和反演效果；

②采集的数据量大是电阻率跨孔CT探测方法的一个重大优势，可以获得关于异常体更多的有效信息，但对大量数据的反演解译处理将非常耗时，亟待提出一种提高反演解译效率的算法；

③不合理的观测方式、不恰当的电极间距和钻孔间距将严重影响电阻率CT方法对孤石目标体的探测和识别能力，导致孤石探测结果存在较大误差甚至探不到孤石的情况发生，因此，亟待提出一种提高孤石探测分辨率和定位精度的科学合理的观测方式；

④在工程实际勘察过程中，既要考虑不同电极间距和孔间距内电阻率CT对孤石的勘探效果，又要考虑用于钻探的较高成本，亟需提出一种较为科学合理同时又较为经济的钻孔平面布置方案，以提高电阻率CT反演成像方法探测孤石的效率。

⑤对于靠近某一钻孔的地质异常体，其反演成像结果中对应的形态会发生较大的畸变，同时会影响临近区域的成像质量，出现“多余构造”或假异常，容易影响或误导地质解释和推断工作。该问题长期困扰三维电阻率跨孔CT探测，目前人们除了依靠经验进行识别和校正之外，没有更有效的方法。影响钻孔附近成像质量的原因主要有以下两点：

1、反演方程对电极近处的异常体的敏感度较高，从而使得电极附近的异常体在反演过程中处于“优势地位”，导致对应的成像结果的形态发生较大畸变；

2、异常体过于靠近电极容易导致电极附近局部电场发生畸变，使成像结果出现偏差。

为此，需要研究并提出一种基于电阻率跨孔CT的地铁盾构区间孤石精细化探测方法，在现有水平和基础上对钻孔中探测电极附近的成像效果进行改进，进一步提高三维电阻率CT的探测和定位精度，为实现地铁盾构掘进机施工隧道条件下隧道掘进面前方孤石体、孤石群的精细化探测提供一条可行的途径。

发明内容

本发明的目的是为克服上述现有技术的不足，提供一种在地铁盾构掘进机施工条件下实现隧道掘进面前方孤石三维空间定位和形态识别、基于电阻率跨孔CT技术的精细化探测方法。该方法采用电阻率跨孔CT方法采集数据，对探测结果进行高精度快速反演解译，实现对孤石体、孤石群的三维空间定位和形态识别。

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

一种基于电阻率跨孔CT的地铁盾构区间孤石精细化探测方法，它包括以下步骤：

(1)根据现场环境特点和工程造价在掘进面前方未开挖段布置钻孔平面；

(2)铺设电阻率跨孔CT测线网，采用二维电阻率跨孔CT探测普查、三维电阻率跨孔CT探测详查的探测方案对地铁隧道盾构掘进前方未开挖段的孤石体进行三维空间定位和形态识别；

(3)将探测电极安装到钻孔中，利用组合式观测模式获取孔间介质的地电结构信息，将采集到的数据提取出来形成一个新的序列，所述组合式观测模式包括四极法电极阵列、三极法电极阵列和二极法电极阵列方法；

(4)对该序列进行综合反演解译，构造携带距离加权函数的三维电阻率反演目标函数及反演方程，并对反演方程进行求解；

(5)形成测区范围内孤石异常体探测的电阻率成像结果图，并结合已有地质分析推断孤石体的大小、空间分布情况。

所述步骤(1)中布置钻孔平面的具体方法为：

沿隧道开挖轴线的两侧布置两排钻孔，钻孔垂直地面打入，将钻孔布置在隧道边界之外，相邻两钻孔间距控制在设定距离以内。

所述步骤(2)的具体方法为：

首先采用二维电阻率跨孔CT探测方法对地铁隧道盾构掘进前方未开挖段的孤石体赋存情况进行普查，圈定疑似单个大直径孤石、孤石群地层或存在安全隐患的重点探查区域，然后在重点探查区域进一步开展三维电阻率跨孔CT探测方法详查，实现孤石体的三维空间定位和形态识别。

所述二维电阻率跨孔CT探测方法为在一对钻孔内布置测线进行数据采集和分析，钻孔连线与隧道轴线垂直或相交，形成孔间二维地电断面成像结果图，能够反映探测剖面内部及外部一定范围内的地质异常；

所述三维电阻率跨孔CT探测方法为在沿隧道轴线两侧各2个相邻钻孔内布置测线进行数据采集和分析，形成孔间三维全空间电阻率成像结果图，能够对钻孔所围成区域内部的地质异常进行精确的三维空间定位和形态识别。

所述步骤(4)中构造的目标函数具体为：

Φ_l＝(Δd-AΔm)^T(Δd-AΔm)+λ(W_lCΔm)^T(W_lCΔm)；

其中，Δd表示实测数据与正演理论数据的差向量，Δm表示模型参数增量向量，A表示敏感度矩阵，C表示光滑度矩阵；λ为拉格朗日常数，用于控制数据方差项和模型方差项之间的权重；w_l表示网格模型中网格的距离加权因子。

所述模型中网格的距离加权因子w_l的计算方法为：

建立直角坐标系，假设四个钻孔P1、P2、P3和P4等间距布置，且距离为a，模型中任一网格i(x,y,z)中心点到各钻孔P1、P2、P3、P4的距离l_i1、l_i2、l_i3、l_i4可表示如下：

$\left\{\begin{array}{c}{l}_{i1}=\sqrt{x^2+y^2}\\ l_{i2}=\sqrt{x^2+{(a-y)}^2}\\ l_{i3}=\sqrt{{(a-x)}^2+{(a-y)}^2}\\ l_{i4}=\sqrt{{(a-x)}^2+y^2}\end{array}\right.;$

则

$w_{\mathrm{li}}(x,y)=1g[\mathrm{\βa}\·(\frac{1}{l_{i1}}+\frac{1}{l_{i2}}+\frac{1}{l_{i3}}+\frac{1}{l_{i4}})];$

上式中，w_li(x,y)表示模型中第i个网格的距离加权因子；a为钻孔间距，β为经验估计常数，通过调节β值可以控制距离加权函数的值域在一定范围内变化；lg表示以10为底的常用对数；

设模型中网格数量为M，将距离加权函数写成矩阵形式，即

其中，w_l1,...,w_lM表示第1～M个网格的距离加权因子。

所述步骤(4)中构造携带距离加权函数的反演方程为：

[A^TA+λ(W_lC)^T(W_lC)]Δm＝A^TΔd

其中，A为偏导数矩阵，Δm为模型参数增量向量，Δd为观测数据，C表示光滑度矩阵；λ为拉格朗日常数，w_l表示网格模型中网格的距离加权因子。

所述偏导数矩阵A为：

$A_{\mathrm{ij}}=\frac{\∂{\mathrm{\ρ}}_{\mathrm{si}}}{\∂{\mathrm{\ρ}}_j}=\underset{\mathrm{\Δ}{\mathrm{\ρ}}_j\→0}{\lim }\frac{\mathrm{\Δ}{\mathrm{\ρ}}_{\mathrm{si}}}{\mathrm{\Δ}{\mathrm{\ρ}}_j}$

其中，A_ij为偏导数矩阵A中第i行第j列个元素，ρ_si为第i个视电阻率值，ρ_j为第j个模型网格的电阻率值，令Δρ_j＝ρ_j-0.99ρ_j，Δρ_si＝ρ_si1-ρ_si2，其中ρ_si1，ρ_si2分别为模型网格电阻率值ρ_j改变前后的视电阻率值。

在对数据进行反演解译分析时采用主从并行模式对偏导数矩阵进行并行解析求解；所述主从并行模式是指在求解偏导数矩阵时采用并行计算，而反演的其他部分依然采用串行计算的算法；主从式并行改造之后，采用分块并行策略，将计算任务平均分配给计算机线程，提高反演计算精度和计算效率。

所述采用主从并行模式对偏导数矩阵进行并行解析求解的具体方法为：

1)将电阻率有限元模型作为参数化模型，模型参数数目为m，设定模型参数初值；

2)在电阻率CT反演中每迭代一次均需要针对每个模型即模型网格电阻率值ρ_j进行正演计算，求取对应模型的地球物理响应；

3)收敛判断：收敛的判据为rus＜ε_inv，所述rus表达式如下：

$\mathrm{rus}=\sqrt{\frac{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\left(\frac{d_{\mathrm{obsi}}-d_{\mathrm{mi}}}{d_{\mathrm{obsi}}}\right)}^2}{N}}$

其中rus为观测数据d_obs与正演理论值d_m之差与观测数据d_obs比值的均方差，d_obsi为第i个实际观测视电阻率数据，d_mi为第i个正演理论观测数据，N为观测数据的个数，ε_inv为反演收敛的容许值；

若rus不满足收敛判据，则转至步骤4)；若rus满足收敛判据，将此时得到的模型参数作为反演的结果输出；

4)偏导数矩阵的并行计算：应用基于共享内存的多处理器多线程并行系统，采用分块并行策略，将每个模型的正演计算任务平均分配给计算机线程；

5)求解反演方程得到新的模型参数，重复步骤2)-步骤3)。

本发明有益效果：

1.本发明提出了一种电阻率跨孔CT二维探测普查、三维探测详查相结合的探测方案，通过普查圈定疑似单个大直径孤石、孤石群地层或存在安全隐患的重点探查区域，然后在重点探查区域详查实现孤石体的三维空间定位和形态识别，大大提高了探测效率和探测精度。

2.本发明提出了一种新型组合式观测模式，综合了常规四极法(AM-BN)、三极法(AM-N)和二极法(A-M)的探测特点和优势，可获得更多有效的地质信息，提高对孤石异常体的探测识别能力。

3.本发明提出了一种三维电阻率CT“距离加权”反演成像方法，施加距离加权约束后反演成像结果中低阻异常的形态与原模型较为吻合，对异常体边界的识别效果明显改善，成像质量较高，表明距离加权约束对于削弱钻孔附近电场畸变影响、解决反演敏感度较高的问题具有较好的效果。

4.本发明提出了一种偏导数矩阵的并行解析快速求解算法，在对数据进行反演解译分析时采用主从并行模式对偏导数矩阵进行并行解析求解，大大提高了数据反演解译的计算精度和计算效率。

附图说明

图1是本发明孤石探测电阻率跨孔CT钻孔平面布置方案示意图；

图2是本发明孤石探测电阻率跨孔CT立体观测示意图；

图3(a)是本发明电阻率跨孔CT四极法观测模式原理图；

图3(b)是本发明电阻率跨孔CT三极法观测模式原理图；

图3(c)是本发明电阻率跨孔CT二极法观测模式原理图；

图4是本发明模型中任一网格与钻孔距离关系示意图；

图5(a)是本发明三维坐标系下距离加权函数分布图；

图5(b)是本发明xy平面上距离加权函数的等值线分布图；

图6是本发明偏导数矩阵并行解析快速求解流程图；

图7(a)是本发明数值算例的孤石体分布模型；

图7(b)是本发明反演得到的电阻率等值面图。

具体实施方式

下面通过具体实例和附图对本发明进行进一步的阐述，应该说明的是，下述说明仅是为了解释本发明，并不对其内容进行限定。

一种基于电阻率跨孔CT的地铁盾构区间孤石精细化探测方法，它包括以下步骤：

A.根据现场环境特点和工程造价在地铁盾构掘进机开挖隧道掘进面前方未开挖段设计科学合理的钻孔平面布置方案，具体如图1所示，沿隧道开挖轴线的两侧布置两排钻孔，钻孔垂直地面打入，尽可能把钻孔布置在隧道边界之外，相邻两钻孔间距一般控制在20m以内。

B.铺设电阻率跨孔CT测线网，采用二维探测普查、三维探测详查相结合的探测方案。首先采用二维电阻率跨孔CT方法对掘进前方未开挖段的孤石体赋存情况进行普查，二维探测是指在一对钻孔内布置测线进行数据采集和分析的探测方法，钻孔连线应与隧道轴线垂直或相交，如图1中的钻孔对P1-P2、P2-P3、P3-P4、P5-P6、P7-P8、P9-P10等，可形成孔间二维地电断面成像结果图，初步圈定疑似单个大直径孤石、孤石群地层或存在安全隐患的重点探查区域(如图1所示，假定在P5～P10围成的区域内发现疑似孤石群地层)；二维探测方法能够较好地揭露孤石异常体的埋藏深度，但无法实现孤石体的三维空间定位和形态识别，而施工中必须在明确孤石体赋存的空间位置和大小后才能进行相应处理，因此需要在重点探查区域进一步开展三维电阻率跨孔CT方法详查，三维探测是指在相邻4个钻孔(沿隧道轴线两侧各2个钻孔)内布置测线进行数据采集和分析的探测方法，如图1中钻孔P5-P6-P7-P8、P7-P8-P9-P10所形成的探测区域，或如图2所示，可形成孔间三维全空间电阻率成像结果图，实现探测区域内部孤石体的三维空间定位和形态识别。

需要说明的是，在进行数据处理时，二维探测区域采用矩形网格进行剖分，三维探测区域采用八节点六面体网格进行剖分。

C.数据采集利用新型组合式观测模式，所述新型组合式观测模式综合了常规四极法(AM-BN)、三极法(AM-N)和二极法(A-M)的探测特点和优势，其中A、B代表供电电极，M、N代表测量电极。基于电阻率跨孔CT的四极法电极阵列数据采集方式如图3(a)所示，设每个钻孔中布置n个电极，钻孔1中电极编号从上到下依次为1#、2#、……、n#，钻孔2中电极编号从上到下依次为n+1#、n+2#、……、2n#，测量分两个循环：

第①循环：A、M位于钻孔1中，N、B位于钻孔2中，测量开始，首先供电电极A位于1#电极，测量电极M位于4#电极，测量电极N位于n+1#电极，供电电极B位于n+4#电极；A、M位置保持不变，N、B逐次分别向下移动一个电极，直至N移至2n-3#电极，同时B移至2n#电极，此时可采集得到(n-3)个数据；A、M同时向下移动一个电极并保持位置不变，N、B按照上述规则逐次移动测量，又得到(n-3)个数据；直至A移至n-3#电极，同时M移至n#电极，第①循环共采集得到(n-3)²个数据。

第②循环：A、M位于钻孔2中，N、B位于钻孔1中，又可得到(n-3)²数据。

因此，基于电阻率跨孔CT的四极法电极阵列一次工作可获得的数据总数为2×(n-3)²个。同理，可根据相同规律得到基于电阻率跨孔CT的三极法和二极法电极阵列采集的数据总数，采集方式分别如图3(a)-3(c)所示,这里不再赘述。

图3(a)～图3(c)三种观测模式叠加形成了本发明新型组合式观测模式。

新型组合式观测模式将3种观测模式采集到的数据提取出来形成一个新的序列，再对该序列进行综合反演解译，使得采集的数据量大大增加，可获得孔间更多有效的地质结构信息，提高了对孤石异常体的探测识别能力。

D.针对孔中探测电极附近成像效果差的问题，提出了“距离加权”的思路，根据反演模型网格距离电极的远近施加不同权重，重点调整钻孔附近区域反演的敏感度以改善其成像效果，对常规光滑约束最小二乘反演进行改进。

常规的三维电阻率探测反演的目标函数为：

Φ＝(Δd-AΔm)^T(Δd-AΔm)+λ(CΔm)^T(CΔm)  (1)

上式中，Δd表示实测数据与正演理论数据的差向量，Δm表示模型参数增量向量，A表示敏感度矩阵，C表示光滑度矩阵；前一项为数据方差项，用于控制实测数据与理论数据的差异最小拟合；后一项为模型方差项，用于控制模型相邻网格电阻率差异最小拟合；λ为拉格朗日常数，用于控制数据方差项和模型方差项之间的权重。

对于八节点六面体三维电阻率有限元网格模型而言，式(1)中光滑矩阵C的表达式如下：

假设C_ij为C中的元素，则

从式(2)可以看出，光滑约束使得相邻网格间电阻率值的光滑过渡，从本质上讲是一种全局宽松约束，因此，单纯施加光滑约束并不能改善探测电极附近成像的效果。

为了克服这一问题，引入“距离加权”约束的思想，根据模型网格与探测电极距离的不同赋予不同的权重，建立权重模型，调整探测区域不同空间位置反演网格的敏感度，对反演的网格电阻率进行加权。建立如图4所示坐标系，假设四个钻孔等间距，距离为a。

则模型中任一网格i(x,y,z)中心点到各钻孔P1、P2、P3、P4的距离l_i1、l_i2、l_i3、l_i4可表示如下：

$\left\{\begin{array}{c}{l}_{i1}=\sqrt{x^2+y^2}\\ l_{i2}=\sqrt{x^2+{(a-y)}^2}\\ l_{i3}=\sqrt{{(a-x)}^2+{(a-y)}^2}\\ l_{i4}=\sqrt{{(a-x)}^2+y^2}\end{array}\right.---\left(3\right)$

在式(3)的基础上，提出了如下距离加权函数表达式：

$w_{\mathrm{li}}(x,y)=1g[\mathrm{\βa}\·(\frac{1}{l_{i1}}+\frac{1}{l_{i2}}+\frac{1}{l_{i3}}+\frac{1}{l_{i4}})]---\left(4\right)$

上式中，w_li(x,y)表示模型中第i个网格的距离加权因子；a为钻孔间距，β为经验估计常数，通过调节β值可以控制距离加权函数的值域在一定范围内变化；lg表示以10为底的常用对数。从公式(4)可以看出，距离加权函数仅仅是x，y的函数，与深度z无关。

分别在三维坐标系下和xy平面上绘制出模型中各网格对应的距离加权因子值分布图像(这里取钻孔间距a＝10m，β＝1.67)，得到的距离加权函数分布图如图5(a)和图5(b)所示，图中红色圆圈代表四个钻孔的位置。从距离加权函数的分布规律中可以看出，只有四个钻孔附近小范围内的模型网格电阻率值得到“有效加权”(加权因子w_li(x,y)＞1)，而处于探测中心区域、距钻孔较远的绝大多数网格电阻率值并未受到距离加权函数的明显影响(w_li(x,y)≈1)。

设模型中网格数量为M，将距离加权函数写成矩阵形式，即

携带距离加权函数的三维电阻率反演目标函数可写成如下形式：

Φ_l＝(Δd-AΔm)^T(Δd-AΔm)+λ(W_lCΔm)^T(W_lCΔm)  (6)

对目标函数公式(6)求极小值，可得到对应的反演方程：

[A^TA+λ(W_lC)^T(W_lC)]Δm＝A^TΔd   (7)

公式(7)即携带“距离加权”约束的三维电阻率CT反演方程。从理论上讲，施加了距离加权约束后，将有效改善反演方程对孔中电极附近区域敏感度较高而造成的成像不稳定问题，提高钻孔附近地质异常体的三维定位和识别精度。与此同时，对于地质异常位于探测范围内其他区域的情况，距离加权约束基本不会对固有的探测成像效果带来影响。

E.在电阻率跨孔CT反演的流程中，最为耗时的是偏导数矩阵的计算和反演方程的求解。为了解决偏导数矩阵计算效率低的瓶颈问题，在保证反演精度的前提下提高反演求解效率，提出了偏导数矩阵并行解析快速求解算法，偏导数矩阵A的计算表达式如下：

$A_{\mathrm{ij}}=\frac{\∂{\mathrm{\ρ}}_{\mathrm{si}}}{\∂{\mathrm{\ρ}}_j}=\underset{\mathrm{\Δ}{\mathrm{\ρ}}_j\→0}{\lim }\frac{\mathrm{\Δ}{\mathrm{\ρ}}_{\mathrm{si}}}{\mathrm{\Δ}{\mathrm{\ρ}}_j}$

上式中A_ij为偏导数矩阵A中第i行第j列个元素，ρ_si为第i个视电阻率值，ρ_j为第j个模型网格的电阻率值。解析法是利用定义将偏导数变差分进行直接求解的方法，令Δρ_j＝ρ_j-0.99ρ_j，Δρ_si＝ρ_si1-ρ_si2，其中ρ_si1，ρ_si2分别为模型网格电阻率值ρ_j改变前后的视电阻率值。

偏导数矩阵并行解析快速求解方法及其流程如图6所示，将电阻率有限元模型作为参数化模型，模型参数数目为m，在电阻率CT反演中每迭代一次均需要针对每个模型参数(即模型网格电阻率值ρ_j)进行正演计算，以求取对应模型的地球物理响应。由于不同模型参数的正演计算是彼此独立的，因此利用主从并行模式进行偏导数矩阵的并行解析快速求解具有天然的合理性。主从并行模式指只在求解偏导数时采用并行计算，而反演的其他部分依然采用串行计算。主从式并行改造之后，采用分块并行策略，将计算任务平均分配给计算机线程。以多核CPU单机计算机为实现平台，应用基于共享内存的多处理器多线程并行系统，使程序自动分配线程进行并行化计算。

具体计算流程如下：

1)将电阻率有限元模型作为参数化模型，模型参数数目为m，设定模型参数初值；

2)在电阻率CT反演中每迭代一次均需要针对每个模型即模型网格电阻率值ρ_j进行正演计算，求取对应模型的地球物理响应；

3)收敛判断：收敛的判据为rus＜ε_inv，所述rus表达式如下：

$\mathrm{rus}=\sqrt{\frac{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\left(\frac{d_{\mathrm{obsi}}-d_{\mathrm{mi}}}{d_{\mathrm{obsi}}}\right)}^2}{N}}$

其中rus为观测数据d_obs与正演理论值d_m之差与观测数据d_obs比值的均方差，d_obsi为第i个实际观测视电阻率数据，d_mi为第i个正演理论观测数据，N为观测数据的个数，ε_inv为反演收敛的容许值，

若rus不满足收敛判据，则转至步骤4)；若rus满足收敛判据，将此时得到的模型参数作为反演的结果输出；

4)偏导数矩阵的并行计算：应用基于共享内存的多处理器多线程并行系统，采用分块并行策略，将每个模型的正演计算任务平均分配给计算机线程；

5)求解反演方程得到新的模型参数，重复步骤2)-步骤3)。

F.经上述步骤对采集到的数据进行反演解译处理后，便可得到测区范围内的孤石异常体电阻率反演成像结果图。图7(a)给出了两个孤石体分布的数值算例模型，图7(b)是采用本发明所述方法综合反演解译得到的电阻率等值面图，从图中可以清晰地看出孤石异常体的埋藏深度、形态和规模，与原始模型基本一致。

通过二维探测普查、三维探测详查相结合的探测方案，可以得到整个测区内较为详细的孤石分布电阻率成像结果图，结合已有地质分析进一步推断孤石体的大小、空间分布等情况，能够为实现地铁盾构掘进机施工隧道安全快速掘进提供有力的技术支持和保障。

上述虽然结合附图对本发明的具体实施方式进行了描述，但并非对本发明保护范围的限制，所属领域技术人员应该明白，在本发明的技术方案的基础上，本领域技术人员不需要付出创造性劳动即可做出的各种修改或变形仍在本发明的保护范围以内。

Claims (10)

1.一种基于电阻率跨孔CT的地铁盾构区间孤石精细化探测方法，其特征在于，它包括以下步骤：

(1)根据现场环境特点和工程造价在掘进面前方未开挖段布置钻孔平面；

(2)铺设电阻率跨孔CT测线网，采用二维电阻率跨孔CT探测普查、三维电阻率跨孔CT探测详查的探测方案对地铁隧道盾构掘进前方未开挖段的孤石体进行三维空间定位和形态识别；

(3)将探测电极安装到钻孔中，利用组合式观测模式获取孔间介质的地电结构信息，将采集到的数据提取出来形成一个新的序列，所述组合式观测模式包括四极法电极阵列、三极法电极阵列和二极法电极阵列方法；

(4)对该序列进行综合反演解译，构造携带距离加权函数的三维电阻率反演目标函数及反演方程，并对反演方程进行求解；

(5)形成测区范围内孤石异常体探测的电阻率成像结果图，并结合已有地质分析推断孤石体的大小、空间分布情况。

2.根据权利要求1所述的方法，其特征在于，所述步骤(1)中布置钻孔平面的具体方法为：

沿隧道开挖轴线的两侧布置两排钻孔，钻孔垂直地面打入，将钻孔布置在隧道边界之外，相邻两钻孔间距控制在设定距离以内。

3.根据权利要求1所述的方法，其特征在于，所述步骤(2)的具体方法为：

首先采用二维电阻率跨孔CT探测方法对地铁隧道盾构掘进前方未开挖段的孤石体赋存情况进行普查，圈定疑似单个大直径孤石、孤石群地层或存在安全隐患的重点探查区域，然后在重点探查区域进一步开展三维电阻率跨孔CT探测方法详查，实现孤石体的三维空间定位和形态识别。

4.根据权利要求3所述的方法，其特征在于，所述二维电阻率跨孔CT探测方法为在一对钻孔内布置测线进行数据采集和分析，钻孔连线与隧道轴线垂直或相交，形成孔间二维地电断面成像结果图，能够反映探测剖面内部及外部一定范围内的地质异常；

所述三维电阻率跨孔CT探测方法为在沿隧道轴线两侧各2个相邻钻孔内布置测线进行数据采集和分析，形成孔间三维全空间电阻率成像结果图，能够对钻孔所围成区域内部的地质异常进行精确的三维空间定位和形态识别。

5.根据权利要求1所述的方法，其特征在于，所述步骤(4)中构造的目标函数具体为：

Φ_l＝(Δd-AΔm)^T(Δd-AΔm)+λ(W_lCΔm)^T(W_lCΔm)；

其中，Δd表示实测数据与正演理论数据的差向量，Δm表示模型参数增量向量，A表示敏感度矩阵，C表示光滑度矩阵；λ为拉格朗日常数，用于控制数据方差项和模型方差项之间的权重；w_l表示网格模型中网格的距离加权因子。

6.根据权利要求5所述的方法，其特征在于，所述模型中网格的距离加权因子w_l的计算方法为：

建立直角坐标系，假设四个钻孔P1、P2、P3和P4等间距布置，且距离为a，模型中任一网格i(x,y,z)中心点到各钻孔P1、P2、P3、P4的距离l_i1、l_i2、l_i3、l_i4可表示如下：

$\left\{\begin{array}{c}{l}_{i1}=\sqrt{x^2+y^2}\\ l_{i2}=\sqrt{x^2+{(a-y)}^2}\\ l_{i3}=\sqrt{{(a-x)}^2+{(a-y)}^2}\\ l_{i4}=\sqrt{{(a-x)}^2+y^2}\end{array}\right.;$

则

$w_{\mathrm{li}}(x,y)=\lg [\mathrm{\βa}\·(\frac{1}{l_{i1}}+\frac{1}{l_{i2}}+\frac{1}{l_{i3}}+\frac{1}{l_{i4}})];$

上式中，w_li(x,y)表示模型中第i个网格的距离加权因子；a为钻孔间距，β为经验估计常数，通过调节β值可以控制距离加权函数的值域在一定范围内变化；lg表示以10为底的常用对数；

设模型中网格数量为M，将距离加权函数写成矩阵形式，即

其中，w_l1,...,w_lM表示第1～M个网格的距离加权因子。

7.根据权利要求1所述的方法，其特征在于，所述步骤(4)中构造携带距离加权函数的反演方程为：

[A^TA+λ(W_lC)^T(W_lC)]Δm＝A^TΔd

其中，A为偏导数矩阵，Δm为模型参数增量向量，Δd为观测数据，C表示光滑度矩阵；λ为拉格朗日常数，w_l表示网格模型中网格的距离加权因子。

8.根据权利要求7所述的方法，其特征在于，所述偏导数矩阵A为：

$A_{\mathrm{ij}}=\frac{\∂{\mathrm{\ρ}}_{\mathrm{si}}}{{\∂\mathrm{\ρ}}_j}=\underset{\mathrm{\Δ}{\mathrm{\ρ}}_j\→0}{\lim }\frac{{\mathrm{\Δ\ρ}}_{\mathrm{si}}}{{\mathrm{\Δ\ρ}}_j}$

其中，A_ij为偏导数矩阵A中第i行第j列个元素，ρ_si为第i个视电阻率值，ρ_j为第j个模型网格的电阻率值，令Δρ_j＝ρ_j-0.99ρ_j，Δρ_si＝ρ_si1-ρ_si2，其中ρ_si1，ρ_si2分别为模型网格电阻率值ρ_j改变前后的视电阻率值。

9.根据权利要求1所述的方法，其特征在于，在对数据进行反演解译分析时采用主从并行模式对偏导数矩阵进行并行解析求解；所述主从并行模式是指在求解偏导数矩阵时采用并行计算，而反演的其他部分依然采用串行计算的算法；主从式并行改造之后，采用分块并行策略，将计算任务平均分配给计算机线程，提高反演计算精度和计算效率。

10.根据权利要求9所述的方法，其特征在于，所述采用主从并行模式对偏导数矩阵进行并行解析求解的具体方法为：

1)将电阻率有限元模型作为参数化模型，模型参数数目为m，设定模型参数初值；

2)在电阻率CT反演中每迭代一次均需要针对每个模型即模型网格电阻率值ρ_j进行正演计算，求取对应模型的地球物理响应；

3)收敛判断：收敛的判据为rus＜ε_inv，所述rus表达式如下：

$\mathrm{rus}=\sqrt{\frac{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\left(\frac{d_{\mathrm{obsi}}-d_{\mathrm{mi}}}{d_{\mathrm{obsi}}}\right)}^2}{N}}$

其中rus为观测数据d_obs与正演理论值d_m之差与观测数据d_obs比值的均方差，d_obsi为第i个实际观测视电阻率数据，d_mi为第i个正演理论观测数据，N为观测数据的个数，ε_inv为反演收敛的容许值；

若rus不满足收敛判据，则转至步骤4)；若rus满足收敛判据，将此时得到的模型参数作为反演的结果输出；

4)偏导数矩阵的并行计算：应用基于共享内存的多处理器多线程并行系统，采用分块并行策略，将每个模型的正演计算任务平均分配给计算机线程；

5)求解反演方程得到新的模型参数，重复步骤2)-步骤3)。