CN103592697B - 隧道等地下工程先验形态信息约束的联合超前预报方法 - Google Patents

Abstract

本发明涉及一种隧道等地下工程先验形态信息约束的联合超前预报方法，其主要步骤为：首先对探测区域收集地质资料、钻探资料分析，对于重点高风险突涌水区域，实施陆地声纳法、瞬变电磁法、地质雷达法、激发极化法超前预报数据采集，选取确定性探测结果异常作为已知先验信息；根据获取的异常体区域，构造该区域的方向向量及其权重，组装成先验信息约束对角阵，对激发极化反演进行约束，形成了先验形态信息约束的联合超前预报方法。该方法实现了联合预报数据约束，显著压制了多解性，提高了综合超前预报的准确性与可靠性。实现了掌子面前方突涌水灾害源的三维空间形态预报，尤其可以实现对突涌水灾害源边界的识别。

Description

隧道等地下工程先验形态信息约束的联合超前预报方法

技术领域

本发明涉及一种隧道或者坑道等地下工程中基于先验形态信息约束的联合超前地质预报方法，具体地说是一种涉及在隧道或坑道中采用电磁类方法（瞬变电磁法、地质雷达法）、弹性波类方法（陆地声纳法、TSP方法）与激发极化法等进行联合超前预报，将电磁类方法与弹性波类方法的探测结果作为确定性的已知先验信息，并施加到激发极化反演目标函数中，从而形成一种隧道等地下工程先验形态信息约束的联合超前预报方法，以实现对隧道前方的不良地质体进行空间定位。

背景技术

隧道施工期突涌水灾害是地下工程施工面临的巨大挑战，隧道超前地质预报对隧道安全施工，避免重大事故发生提供了重要的保障，近年来发挥着越来越重要的作用。目前隧道超前预报方法主要采用地球物理勘探方法，由于地球物理勘探存在多解性的固有难题，单一的预报方法往往导致预报结果的偏差，甚至错误。反射类地震波方法（TSP等），均对岩性界面、断层、岩溶空洞等有较好的反映；电磁类预报方法（如地质雷达法、瞬变电磁法等）对水体较为敏感，可以用来预报工作面前方的含水体；激发极化对突涌水含水构造响应敏感，且能估算含水构造的涌水量。为此结合不同方法的特点，很多学者，采用弹性波类方法、电磁类方法与电法等地球物理探测方法，提出了隧道联合超前预报方法。

但目前的综合超前预报仍面临以下问题：①目前的综合地球物理探测大多是在数据解释阶段对多种探测方法的结果进行人为的比对分析，探测可靠性依靠个人经验，未实现多种探测方法数据的相互约束，没有发挥不同探测方法的优势，多解性仍严重。②目前综合超前预报探测突涌水灾害源水体普遍采用电磁法，但由于隧道狭小复杂环境，往往不能对水体三维形态、位置进行准确预报，更无法反映含水体的边界等重要特征信息。

发明内容

本发明的目的是为了克服上述现有技术的不足，提供一种隧道等地下工程先验形态信息约束的联合超前预报方法，它采用多种超前预报结果信息作为先验约束，实现了以激发极化法为载体的先验形态约束联合超前预报方法，显著压制了解译的多解性，可以反映突涌水灾害源的三维形态与边界，提高了联合超前预报的准确性与可靠性。

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

一种隧道等地下工程先验形态信息约束的联合超前预报方法，包括如下步骤：

A.首先收集探测区域的地质资料、钻探资料；

B.对于重点高风险突涌水区域，实施陆地声纳法、瞬变电磁法、地质雷达法、激发极化法进行超前预报数据采集；

C.选取陆地声纳法、瞬变电磁法与地质雷达法探测结果中的异常信息作为已知先验信息；

D.建立激发极化法的三维反演坐标***，并将步骤C中各异常信息对应到激发极化法的三维反演坐标***中，从而获取陆地声纳、瞬变电磁与地质雷达探测的异常体区域信息；

E.利用空间形态先验约束方法，根据获取的异常体区域，构造该区域的方向向量及其权重，组装成先验信息约束对角阵，将已知先验形态信息施加到激发极化法反演方程中；

F.用先验形态信息约束反演方程，进行反演迭代，得到掌子面前方突涌水灾害源的三维空间形态与边界信息，实现突涌水灾害源的联合超前预报。

所述步骤C中，对陆地声纳法、瞬变电磁法与地质雷达法进行解释，选取确定性探测结果异常作为已知先验信息，陆地声纳法对断层、溶洞等界面探测效果较好，将陆地声纳探测的异常界面区域作为已知先验信息，瞬变电磁与地质雷达对低电阻区域探测效果好，将瞬变电磁与地质雷达探测的异常区域作为已知先验信息。

所述步骤D中，建立激发极化法的三维反演坐标***，沿着掌子面宽度方向为X轴、高度方向为Y轴、隧道前方为Z轴，将陆地声纳探测结果的波速异常区域的坐标按隧道实际坐标一一对应到激发极化坐标***，将瞬变电磁探测结果的低电阻率异常区域的坐标亦按照隧道实际坐标一一对应到激发极化坐标***，同理，将地质雷达探测的异常区域仍按照隧道实际坐标一一对应到激发极化坐标***，从而获取陆地声纳、瞬变电磁与地质雷达探测的异常体区域信息。

所述步骤E中，根据获取陆地声纳、瞬变电磁与地质雷达探测的异常体区域，采用三个正交方向向量即法向量及两个切向量及各向异性的权重系数表征已知先验形态信息获取的异常体边界，具体做法是，选坐标轴XYZ方向单位向量为三个方向向量，沿着异常区域边界电阻率变化比较缓慢，该方向的方向向量权重较大，一般取1.0-2.0；跨过异常区域边界的方向电阻率变化较大，该方向的方向向量权重较小，一般取0.01-0.2，从而控制激发极化反演中异常体边界区域的法方向及切方向电阻率差异程度，组装成先验信息约束对角阵，将已知先验形态信息施加到激发极化法反演方程中。

具体做法是，公式（1）为目标函数，对公式（1）进行变分得到反演方程（2）

$\begin{array}{c}\mathrm{\Φ}={(\mathrm{\Δd}-\mathrm{A\Δm})}^T(\mathrm{\Δd}-\mathrm{A\Δm})+\mathrm{\Δ}{m}^T{C}_x^T{\mathrm{\Λ}}_x{C}_x\mathrm{\Δm}+\mathrm{\Δ}{m}^T{C}_y^T{\mathrm{\Λ}}_y{C}_y\mathrm{\Δm}+\mathrm{\Δ}{m}^T{C}_z^T{\mathrm{\Λ}}_z{C}_z\mathrm{\Δm}\\ +\mathrm{\Δ}{m}^T{C}_x^T{\mathrm{\Λ}}_{\mathrm{xy}}{C}_y\mathrm{\Δm}+\mathrm{\Δ}{m}^T{C}_x^T{\mathrm{\Λ}}_{\mathrm{xz}}{C}_z\mathrm{\Δm}+\mathrm{\Δ}{m}^T{C}_y^T{\mathrm{\Λ}}_{\mathrm{yz}}{C}_z\mathrm{\Δm}\end{array}---\left(1\right)$

(A^TA+W^TΛ₁W+P^TΛ₂Q)Δm＝A^TΔd (2)

其中，

$\begin{array}{c}{W}^T{\mathrm{\Λ}}_{1}W=C_x^T{\mathrm{\Λ}}_x{C}_x+C_y^T{\mathrm{\Λ}}_y{C}_y+C_z^T{\mathrm{\Λ}}_z{C}_z\\ P^T{\mathrm{\Λ}}_{2}Q=C_x^T{\mathrm{\Λ}}_{\mathrm{xy}}{C}_y+C_y^T{\mathrm{\Λ}}_{\mathrm{xy}}{C}_x+C_x^T{\mathrm{\Λ}}_{\mathrm{xz}}{C}_z+C_z^T{\mathrm{\Λ}}_{\mathrm{xz}}{C}_x+C_y^T{\mathrm{\Λ}}_{\mathrm{yz}}{C}_z+C_z^T{\mathrm{\Λ}}_{\mathrm{yz}}{C}_y\end{array}---\left(3\right)$

式中Φ为反演目标函数，Δd为观测数据的差值，A为偏导数矩阵，Δm为网格的模型参数修正量，C_x、C_y、C_z分别为X、Y、Z方向的光滑度矩阵，Λ_x、Λ_y、Λ_z、Λ_xy、Λ_xz、Λ_yz为携带着方向向量及其权重的对角阵；

$\begin{array}{c}{\mathrm{\Λ}}_{\mathrm{xi}}=\frac{\left|V_i\right|}{{\left(\mathrm{\Δ}{x}_i\right)}^2}[{\left({\mathrm{\λ}}_1^{\left(i\right)}v{1}_1^{\left(i\right)}\right)}^2+{\left({\mathrm{\λ}}_2^{\left(i\right)}v{1}_2^{\left(i\right)}\right)}^2+{\left({\mathrm{\λ}}_3^{\left(i\right)}v{1}_3^{\left(i\right)}\right)}^2]\\ {\mathrm{\Λ}}_{\mathrm{yi}}=\frac{\left|V_i\right|}{{\left(\mathrm{\Δ}{y}_i\right)}^2}[{\left({\mathrm{\λ}}_1^{\left(i\right)}v{2}_1^{\left(i\right)}\right)}^2+{\left({\mathrm{\λ}}_2^{\left(i\right)}v{2}_2^{\left(i\right)}\right)}^2+{\left({\mathrm{\λ}}_3^{\left(i\right)}v{2}_3^{\left(i\right)}\right)}^2]\\ {\mathrm{\Λ}}_{\mathrm{yi}}=\frac{\left|V_i\right|}{{\left(\mathrm{\Δ}{z}_i\right)}^2}[{\left({\mathrm{\λ}}_1^{\left(i\right)}v{3}_1^{\left(i\right)}\right)}^2+{\left({\mathrm{\λ}}_2^{\left(i\right)}v{3}_2^{\left(i\right)}\right)}^2+{\left({\mathrm{\λ}}_3^{\left(i\right)}v{3}_3^{\left(i\right)}\right)}^2]\\ {\mathrm{\Λ}}_{\mathrm{xyi}}=\frac{\left|V_i\right|}{\mathrm{\Δ}{x}_i\mathrm{\Δ}{y}_i}[{\left({\mathrm{\λ}}_1^{\left(i\right)}\right)}^{2}v{1}_1^{\left(i\right)}v{2}_1^{\left(i\right)}+{\left({\mathrm{\λ}}_2^{\left(i\right)}\right)}^{2}v{1}_2^{\left(i\right)}v{2}_2^{\left(i\right)}+{\left({\mathrm{\λ}}_3^{\left(i\right)}\right)}^{2}v{1}_3^{\left(i\right)}v{2}_3^{\left(i\right)}]\\ {\mathrm{\Λ}}_{\mathrm{xzi}}=\frac{\left|V_i\right|}{\mathrm{\Δ}{x}_i\mathrm{\Δ}{z}_i}[{\left({\mathrm{\λ}}_1^{\left(i\right)}\right)}^{2}v{1}_1^{\left(i\right)}v{3}_1^{\left(i\right)}+{\left({\mathrm{\λ}}_2^{\left(i\right)}\right)}^{2}v{1}_2^{\left(i\right)}v{3}_2^{\left(i\right)}+{\left({\mathrm{\λ}}_3^{\left(i\right)}\right)}^{2}v{1}_3^{\left(i\right)}v{3}_3^{\left(i\right)}]\\ {\mathrm{\Λ}}_{\mathrm{yzi}}=\frac{\left|V_i\right|}{\mathrm{\Δ}{y}_i\mathrm{\Δ}{z}_i}[{\left({\mathrm{\λ}}_1^{\left(i\right)}\right)}^{2}v{2}_1^{\left(i\right)}v{3}_1^{\left(i\right)}+{\left({\mathrm{\λ}}_2^{\left(i\right)}\right)}^{2}v{2}_2^{\left(i\right)}v{3}_2^{\left(i\right)}+{\left({\mathrm{\λ}}_3^{\left(i\right)}\right)}^{2}v{2}_3^{\left(i\right)}v{3}_3^{\left(i\right)}]\end{array}---\left(4\right)$

为第i个单元的k（k=1，2，3）个向量的元素， 为第i个单元的三个方向向量权重，V_i为第i个单元的体积，Δx_i、Δy_i、Δz_i为第i个单元沿x、y、z方向的长度。

本发明的有益效果是：本发明提出了一种隧道等地下工程先验形态信息约束的联合超前预报方法，通过方向向量权重组成的空间形态对角阵，将地质雷达法、瞬变电磁法、陆地声纳法探测结果作为已知先验信息，形成了以激发极化法为载体的先验形态约束联合超前预报方法，实现了掌子面前方突涌水灾害源的三维空间形态预报，尤其可以实现对突涌水灾害源边界的识别，具有以下特色：

1>提出了采用地质雷达法、瞬变电磁法、陆地声纳法等预报方法获取的已知先验形态信息作为激发极化反演约束的思想，实现了多种超前预报方法的数据融合约束，解决了传统联合超前预报简单组合未实现数据之间约束的问题。

2>通过三个正交方向向量(法向量及两个切向量)及各向异性的权重系数表征已知先验形态信息边界，形成空间形态约束的联合超前预报方法，提高了联合超前预报的准确性与可靠性。

3>通过方向向量权重组成的空间形态对角阵，分别控制X、Y、Z方向的光滑度矩阵C_x、C_y、C_z，将其他预报方法获取的空间形态等先验信息施加到激发极化反演方程中，解决了掌子面前方突涌水灾害源的三维空间形态预报与边界识别的关键难题。

附图说明

图1是本发明联合超前预报方法流程图。

图2是空间先验形态约束原型

图3是本发明空间形态约束联合预报方法反演结果。

图4是本发明空间形态约束联合预报工程应用结果。

具体实施方式

下面结合附图与实施例对本发明做进一步说明。

图1中，本发明的先验形态信息约束的联合超前预报方法，将原来各种预报方法简单比对分析的联合预报，改变为采用其他预报方法确定性的结果约束激发极化法反演的联合超前预报方法，本发明采集数据的具体流程是：

A.首先对探测区域收集地质资料、钻探资料，进行地质分析推断初步推断是否存在突涌水灾害构造，如果存在突涌水灾害构造，则选出重点高风险突涌水区域；

B.对于重点高风险突涌水区域，实施陆地声纳法、瞬变电磁法、地质雷达法、激发极化法超前预报数据采集；

C.对陆地声纳法、瞬变电磁法与地质雷达法进行解释，选取确定性探测结果异常作为已知先验信息，陆地声纳法对断层、溶洞等界面探测效果较好，将陆地声纳探测的异常界面区域作为已知先验信息，瞬变电磁与地质雷达对低电阻区域探测效果好，将瞬变电磁与地质雷达探测的异常区域作为已知先验信息；

D.建立激发极化法的三维反演坐标***，沿着掌子面宽度方向为X轴、高度方向为Y轴、隧道前方为Z轴，将陆地声纳探测结果的波速异常区域的坐标按隧道实际坐标一一对应到激发极化坐标***，将瞬变电磁探测结果的低电阻率异常区域的坐标亦按照隧道实际坐标一一对应到激发极化坐标***，同理，将地质雷达探测的异常区域仍按照隧道实际坐标一一对应到激发极化坐标***，从而获取陆地声纳、瞬变电磁与地质雷达探测的异常体区域信息；

E.根据获取陆地声纳、瞬变电磁与地质雷达探测的异常体区域，采用三个正交方向向量即法向量及两个切向量及各向异性的权重系数表征已知先验形态信息获取的异常体边界，具体做法是，选坐标轴XYZ方向单位向量为三个方向向量，沿着异常区域边界电阻率变化比较缓慢，该方向的方向向量权重较大，一般取1.0-2.0；跨过异常区域边界的方向电阻率变化较大，该方向的方向向量权重较小，一般取0.01-0.2，从而控制激发极化反演中异常体边界区域的法方向及切方向电阻率差异程度，按照公式（4）组装成先验信息约束对角阵，将已知先验形态信息施加到激发极化法反演方程中，如下式；具体做法是，公式（1）为目标函数，对公式（1）进行变分可得到反演方程（2）。

$\begin{array}{c}\mathrm{\Φ}={(\mathrm{\Δd}-\mathrm{A\Δm})}^T(\mathrm{\Δd}-\mathrm{A\Δm})+\mathrm{\Δ}{m}^T{C}_x^T{\mathrm{\Λ}}_x{C}_x\mathrm{\Δm}+\mathrm{\Δ}{m}^T{C}_y^T{\mathrm{\Λ}}_y{C}_y\mathrm{\Δm}+\mathrm{\Δ}{m}^T{C}_z^T{\mathrm{\Λ}}_z{C}_z\mathrm{\Δm}\\ +\mathrm{\Δ}{m}^T{C}_x^T{\mathrm{\Λ}}_{\mathrm{xy}}{C}_y\mathrm{\Δm}+\mathrm{\Δ}{m}^T{C}_x^T{\mathrm{\Λ}}_{\mathrm{xz}}{C}_z\mathrm{\Δm}+\mathrm{\Δ}{m}^T{C}_y^T{\mathrm{\Λ}}_{\mathrm{yz}}{C}_z\mathrm{\Δm}\end{array}---\left(1\right)$

(A^TA+W^TΛ₁W+P^TΛ₂Q)Δm＝A^TΔd (2)

其中，

$\begin{array}{c}{W}^T{\mathrm{\Λ}}_{1}W=C_x^T{\mathrm{\Λ}}_x{C}_x+C_y^T{\mathrm{\Λ}}_y{C}_y+C_z^T{\mathrm{\Λ}}_z{C}_z\\ P^T{\mathrm{\Λ}}_{2}Q=C_x^T{\mathrm{\Λ}}_{\mathrm{xy}}{C}_y+C_y^T{\mathrm{\Λ}}_{\mathrm{xy}}{C}_x+C_x^T{\mathrm{\Λ}}_{\mathrm{xz}}{C}_z+C_z^T{\mathrm{\Λ}}_{\mathrm{xz}}{C}_x+C_y^T{\mathrm{\Λ}}_{\mathrm{yz}}{C}_z+C_z^T{\mathrm{\Λ}}_{\mathrm{yz}}{C}_y\end{array}---\left(3\right)$

式中Φ为反演目标函数，Δd为观测数据的差值，A为偏导数矩阵，Δm为网格的模型参数修正量，C_x、C_y、C_z分别为X、Y、Z方向的光滑度矩阵，Λ_x、Λ_y、Λ_z、Λ_xy、Λ_xz、Λ_yz为携带着方向向量及其权重的对角阵。

$\begin{array}{c}{\mathrm{\Λ}}_{\mathrm{xi}}=\frac{\left|V_i\right|}{{\left(\mathrm{\Δ}{x}_i\right)}^2}[{\left({\mathrm{\λ}}_1^{\left(i\right)}v{1}_1^{\left(i\right)}\right)}^2+{\left({\mathrm{\λ}}_2^{\left(i\right)}v{1}_2^{\left(i\right)}\right)}^2+{\left({\mathrm{\λ}}_3^{\left(i\right)}v{1}_3^{\left(i\right)}\right)}^2]\\ {\mathrm{\Λ}}_{\mathrm{yi}}=\frac{\left|V_i\right|}{{\left(\mathrm{\Δ}{y}_i\right)}^2}[{\left({\mathrm{\λ}}_1^{\left(i\right)}v{2}_1^{\left(i\right)}\right)}^2+{\left({\mathrm{\λ}}_2^{\left(i\right)}v{2}_2^{\left(i\right)}\right)}^2+{\left({\mathrm{\λ}}_3^{\left(i\right)}v{2}_3^{\left(i\right)}\right)}^2]\\ {\mathrm{\Λ}}_{\mathrm{yi}}=\frac{\left|V_i\right|}{{\left(\mathrm{\Δ}{z}_i\right)}^2}[{\left({\mathrm{\λ}}_1^{\left(i\right)}v{3}_1^{\left(i\right)}\right)}^2+{\left({\mathrm{\λ}}_2^{\left(i\right)}v{3}_2^{\left(i\right)}\right)}^2+{\left({\mathrm{\λ}}_3^{\left(i\right)}v{3}_3^{\left(i\right)}\right)}^2]\\ {\mathrm{\Λ}}_{\mathrm{xyi}}=\frac{\left|V_i\right|}{\mathrm{\Δ}{x}_i\mathrm{\Δ}{y}_i}[{\left({\mathrm{\λ}}_1^{\left(i\right)}\right)}^{2}v{1}_1^{\left(i\right)}v{2}_1^{\left(i\right)}+{\left({\mathrm{\λ}}_2^{\left(i\right)}\right)}^{2}v{1}_2^{\left(i\right)}v{2}_2^{\left(i\right)}+{\left({\mathrm{\λ}}_3^{\left(i\right)}\right)}^{2}v{1}_3^{\left(i\right)}v{2}_3^{\left(i\right)}]\\ {\mathrm{\Λ}}_{\mathrm{xzi}}=\frac{\left|V_i\right|}{\mathrm{\Δ}{x}_i\mathrm{\Δ}{z}_i}[{\left({\mathrm{\λ}}_1^{\left(i\right)}\right)}^{2}v{1}_1^{\left(i\right)}v{3}_1^{\left(i\right)}+{\left({\mathrm{\λ}}_2^{\left(i\right)}\right)}^{2}v{1}_2^{\left(i\right)}v{3}_2^{\left(i\right)}+{\left({\mathrm{\λ}}_3^{\left(i\right)}\right)}^{2}v{1}_3^{\left(i\right)}v{3}_3^{\left(i\right)}]\\ {\mathrm{\Λ}}_{\mathrm{yzi}}=\frac{\left|V_i\right|}{\mathrm{\Δ}{y}_i\mathrm{\Δ}{z}_i}[{\left({\mathrm{\λ}}_1^{\left(i\right)}\right)}^{2}v{2}_1^{\left(i\right)}v{3}_1^{\left(i\right)}+{\left({\mathrm{\λ}}_2^{\left(i\right)}\right)}^{2}v{2}_2^{\left(i\right)}v{3}_2^{\left(i\right)}+{\left({\mathrm{\λ}}_3^{\left(i\right)}\right)}^{2}v{2}_3^{\left(i\right)}v{3}_3^{\left(i\right)}]\end{array}---\left(4\right)$

为第i个单元的k（k=1，2，3）个向量的元素， 为第i个单元的三个方向向量权重，V_i为第i个单元的体积，Δx_i、Δy_i、Δz_i为第i个单元沿x、y、z方向的长度。

F.对先验形态信息约束反演方程（2），进行反演迭代，对携带空间形态先验信息约束的激发极化方程进行反演，得到掌子面前方突涌水灾害源的三维空间形态与边界信息，实现突涌水灾害源的联合超前预报。

图2、图3中，本发明空间形态先验约束的算例，在地面下存在两个倾斜异常构造与一个方形异常构造，设定正常区域电阻率为1000Ω.m，异常区域电阻率为10Ω.m，如图2所示。假设通过地质雷达等其他地球物理预报方法获得了异常体的空间形态信息，在本例中空间形态信息为两个倾斜构造与方形构造的边界界面区域，三个方向向量选择沿坐标轴XYZ方向的单位向量（1，0，0）、（0，1，0）与（0，0，1），沿着异常体边界的权重选择0.02，跨越异常体边界的权重选择1.0，组装成先验信息约束对角阵公式（4），并将空间形态信息施加到式（1）中，经过反演迭代，反演结果如图3所示。图3中，携带空间形态先验约束的联合预报可以较准确的反映异常体的位置，未出现多余构造，异常体与正常区域的边界较为清晰，便于分辨和识别，且可以反映异常体的空间走向。

图4中，本发明在某隧道工程超前探测的工程案例，在此次联合超前预报中采用了地质雷达法、地震波反射类方法、激发极化法。地震波反射类方法与地质雷达法确定的异常体边界，但不能确定异常体是否含水是否为突涌水灾害源。为此，为实现多种预报方法数据的约束融合，将确定的已知先验形态信息施加到激发极化反演中，可以得到边界清晰的导水裂隙，实现了突涌水灾害源的三维空间形态定位与边界识别。

Claims (5)

1.一种隧道或坑道地下工程先验形态信息约束的联合超前预报方法，其特征是，包括如下步骤：

A.首先收集探测区域的地质资料、钻探资料；

B.对于重点高风险突涌水区域，实施陆地声纳法、瞬变电磁法、地质雷达法、激发极化法进行超前预报数据采集；

C.选取陆地声纳法、瞬变电磁法与地质雷达法探测结果中的异常信息作为已知先验信息；

D.建立激发极化法的三维反演坐标***，并将步骤C中各异常信息对应到激发极化法的三维反演坐标***中，从而获取陆地声纳、瞬变电磁与地质雷达探测的异常体区域信息；

E.利用空间形态先验约束方法，根据获取的异常体区域，构造该区域的方向向量及其权重，组装成先验信息约束对角阵，将已知先验形态信息施加到激发极化法反演方程中；

F.用先验形态信息约束反演方程，进行反演迭代，得到掌子面前方突涌水灾害源的三维空间形态与边界信息，实现突涌水灾害源的联合超前预报。

2.如权利要求1所述的隧道或坑道地下工程先验形态信息约束的联合超前预报方法，其特征是，所述步骤C中，对陆地声纳法、瞬变电磁法与地质雷达法进行解释，选取确定性探测结果异常作为已知先验信息，陆地声纳法对断层、溶洞的界面探测效果较好，将陆地声纳探测的异常界面区域作为已知先验信息，瞬变电磁与地质雷达对低电阻区域探测效果好，将瞬变电磁与地质雷达探测的异常区域作为已知先验信息。

3.如权利要求1所述的隧道或坑道地下工程先验形态信息约束的联合超前预报方法，其特征是，所述步骤D中，建立激发极化法的三维反演坐标***，沿着掌子面宽度方向为X轴、高度方向为Y轴、隧道前方为Z轴，将陆地声纳探测结果的波速异常区域的坐标按隧道实际坐标一一对应到激发极化坐标***，将瞬变电磁探测结果的低电阻率异常区域的坐标亦按照隧道实际坐标一一对应到激发极化坐标***，同理，将地质雷达探测的异常区域仍按照隧道实际坐标一一对应到激发极化坐标***，从而获取陆地声纳、瞬变电磁与地质雷达探测的异常体区域信息。

4.如权利要求1所述的隧道或坑道地下工程先验形态信息约束的联合超前预报方法，其特征是，所述步骤E中，根据获取陆地声纳、瞬变电磁与地质雷达探测的异常体区域，采用三个正交方向向量即法向量及两个切向量及各向异性的权重系数表征已知先验形态信息获取的异常体边界，具体做法是，选坐标轴XYZ方向单位向量为三个方向向量，沿着异常区域边界电阻率变化比较缓慢，该方向的方向向量权重较大，取1.0-2.0；跨过异常区域边界的方向电阻率变化较大，该方向的方向向量权重较小，取0.01-0.2，从而控制激发极化反演中异常体边界区域的法方向及切方向电阻率差异程度，组装成先验信息约束对角阵，将已知先验形态信息施加到激发极化法反演方程中。

5.如权利要求4所述的隧道或坑道地下工程先验形态信息约束的联合超前预报方法，其特征是，具体做法是，公式(1)为目标函数，对公式(1)进行变分得到反演方程(2)

$\begin{array}{c}\mathrm{\Φ}={(\mathrm{\Δ}d-A\mathrm{\Δ}m)}^T(\mathrm{\Δ}d-A\mathrm{\Δ}m)+{\mathrm{\Δm}}^T{C}_x^T{\mathrm{\Λ}}_x{C}_x\mathrm{\Δ}m+{\mathrm{\Δm}}^T{C}_y^T{\mathrm{\Λ}}_y{C}_y\mathrm{\Δ}m+{\mathrm{\Δm}}^T{C}_z^T{\mathrm{\Λ}}_z{C}_z\mathrm{\Δ}m\\ +{\mathrm{\Δm}}^T{C}_x^T{\mathrm{\Λ}}_{xy}{C}_y\mathrm{\Δ}m+{\mathrm{\Δm}}^T{C}_x^T{\mathrm{\Λ}}_{xz}{C}_z\mathrm{\Δ}m+{\mathrm{\Δm}}^T{C}_y^T{\mathrm{\Λ}}_{yz}{C}_z\mathrm{\Δ}m\end{array}---\left(1\right)$

(A^TA+W^TΛ₁W+P^TΛ₂Q)Δm＝A^TΔd (2)

其中，

$\begin{array}{c}{W}^T{\mathrm{\Λ}}_{1}W=C_x^T{\mathrm{\Λ}}_x{C}_x+C_y^T{\mathrm{\Λ}}_y{C}_y+C_z^T{\mathrm{\Λ}}_z{C}_z\\ P^T{\mathrm{\Λ}}_{2}Q=C_x^T{\mathrm{\Λ}}_{xy}{C}_y+C_y^T{\mathrm{\Λ}}_{xy}{C}_x+C_x^T{\mathrm{\Λ}}_{xy}{C}_z+C_z^T{\mathrm{\Λ}}_{xz}{C}_x+C_y^T{\mathrm{\Λ}}_{yz}{C}_z+C_z^T{\mathrm{\Λ}}_{yz}{C}_y\end{array}---\left(3\right)$

式中Φ为反演目标函数，Δd为观测数据的差值，A为偏导数矩阵，Δm为网格的模型参数修正量，C_x、C_y、C_z分别为X、Y、Z方向的光滑度矩阵，Λ_x、Λ_y、Λ_z、Λ_xy、Λ_xz、Λ_yz为携带着方向向量及其权重的对角阵；

$\begin{array}{c}{\mathrm{\Λ}}_{xi}=\frac{\left|V_i\right|}{{\left({\mathrm{\Δx}}_i\right)}^2}\[{\left({\mathrm{\λ}}_1^{\left(i\right)}v{1}_1^{\left(i\right)}\right)}^2+{\left({\mathrm{\λ}}_2^{\left(i\right)}v{1}_2^{\left(i\right)}\right)}^2+{\left({\mathrm{\λ}}_3^{\left(i\right)}v{1}_3^{\left(i\right)}\right)}^2\]\\ {\mathrm{\Λ}}_{yi}=\frac{\left|V_i\right|}{{\left({\mathrm{\Δy}}_i\right)}^2}\[{\left({\mathrm{\λ}}_1^{\left(i\right)}v{2}_1^{\left(i\right)}\right)}^2+{\left({\mathrm{\λ}}_2^{\left(i\right)}v{2}_2^{\left(i\right)}\right)}^2+{\left({\mathrm{\λ}}_3^{\left(i\right)}v{2}_3^{\left(i\right)}\right)}^2\]\\ {\mathrm{\Λ}}_{yi}=\frac{\left|V_i\right|}{{\left({\mathrm{\Δz}}_i\right)}^2}\[{\left({\mathrm{\λ}}_1^{\left(i\right)}v{3}_1^{\left(i\right)}\right)}^2+{\left({\mathrm{\λ}}_2^{\left(i\right)}v{3}_2^{\left(i\right)}\right)}^2+{\left({\mathrm{\λ}}_3^{\left(i\right)}v{3}_3^{\left(i\right)}\right)}^2\]\\ {\mathrm{\Λ}}_{xyi}=\frac{\left|V_i\right|}{{\mathrm{\Δx}}_i{\mathrm{\Δy}}_i}\[{\left({\mathrm{\λ}}_1^{\left(i\right)}\right)}^{2}v{1}_1^{\left(i\right)}v{2}_1^{\left(i\right)}+{\left({\mathrm{\λ}}_2^{\left(i\right)}\right)}^{2}v{1}_2^{\left(i\right)}v{2}_2^{\left(i\right)}+{\left({\mathrm{\λ}}_3^{\left(i\right)}\right)}^{2}v{1}_3^{\left(i\right)}v{2}_3^{\left(i\right)}\]\\ {\mathrm{\Λ}}_{xyi}=\frac{\left|V_i\right|}{{\mathrm{\Δx}}_i{\mathrm{\Δz}}_i}\[{\left({\mathrm{\λ}}_1^{\left(i\right)}\right)}^{2}v{1}_1^{\left(i\right)}v{3}_1^{\left(i\right)}+{\left({\mathrm{\λ}}_2^{\left(i\right)}\right)}^{2}v{1}_2^{\left(i\right)}v{3}_2^{\left(i\right)}+{\left({\mathrm{\λ}}_3^{\left(i\right)}\right)}^{2}v{1}_3^{\left(i\right)}v{3}_3^{\left(i\right)}\]\\ {\mathrm{\Λ}}_{yzi}=\frac{\left|V_i\right|}{{\mathrm{\Δy}}_i{\mathrm{\Δz}}_i}\[{\left({\mathrm{\λ}}_1^{\left(i\right)}\right)}^{2}v{2}_1^{\left(i\right)}v{3}_1^{\left(i\right)}+{\left({\mathrm{\λ}}_2^{\left(i\right)}\right)}^{2}v{2}_2^{\left(i\right)}v{3}_2^{\left(i\right)}+{\left({\mathrm{\λ}}_3^{\left(i\right)}\right)}^{2}v{2}_3^{\left(i\right)}v{3}_3^{\left(i\right)}\]\end{array}---\left(4\right)$

为第i个单元的k个向量的元素，k＝1，2，3； 为第i个单元的三个方向向量权重，V_i为第i个单元的体积，Δx_i、Δy_i、Δz_i为第i个单元沿x、y、z方向的长度。