基于多分量观测系统的城市工程地震探测综合成像方法与流程

本发明属于工程地震数据成像处理领域，具体涉及一种城市工程地震探测的成像方法。

背景技术：

随着我国国民经济建设的飞速发展，城市化进程不断深入，城市工程建设在规划、设计、施工阶段都必须对建设区域内的地质情况及地下埋设物有一个系统的勘察，在建设工程中及建成后还必须对工程质量进行检测和监测。工程物探是用于探查工程地质问题(如岩土分层、基岩面形态、岩溶分布、断层裂隙、孤石及复杂地基土等)、环境地质问题(如地下采空塌陷、地面沉降、废弃箱涵、暗河与管线等)的有效手段，其在城市工程建设的作用越来越重要，特别是工程地震勘探技术，由于其采用弹性波方法，具有抗金属与电器干扰能力强，探测深度较大的显著特征，是目前城市工程物探中急需发展的方法。

工程地震勘探在城市地下探测的应用，对工程地震数据处理的精度提出了更高的要求，其数据处理结果将直接影响后期资料解释的准确性和可靠度。目前，工程地震勘探的方法主要是建立在野外石油、煤炭等能源地震勘探的基础上，其数据处理一般包括前期的预处理、数字滤波、速度分析、叠加偏移等，特别是地震偏移技术，它是用计算机按一定的计算方法对观测数据进行处理，使之成为反映地下地质分层界面位置及反射系数值的反射界面的像，是地震勘探数据处理的关键技术，其在一定程度上代表了地震勘探技术的发展程度。绕射扫描偏移是建立在射线偏移基础上使反射波自动归位到真实位置上的一种方法，根据惠更斯原理，地下每一个反射点G都可以看成是一个子波震源。进行绕射扫描偏移时，对计算区域进行网格剖分，把每一个网格点看成是一个反射点，当对探测平面X-Z按一定网格大小划分的每一点G都进行计算，只要划分的足够细，总可以在所要求的精度上反映反射点的全部可能位置。这样，使反射界面上的叠加扫描点G的总振幅相对增大，不在反射界面上的扫描点G的总振幅相对减小，既提高了信噪比，又把反射界面自动偏移到其空间真实位置上去。

城市工程地震勘探有其自身特点，其探测的主要对象位于地下几米到几十米范围内，环境噪声及地面车辆干扰严重，区域范围地形复杂，受人类社会活动影响较大，地下结构横向变化较多，探测区间面积较小，观测系统布置较小，覆盖次数不足，但同时它对探测效率及探测精度要求更高，需要在有限的空间范围内快速高效的完成数据采集，对成像剖面结果要求更高的分辨率，更能反映真实的地下结构及异常位置，因此，工程地震数据需要更严格更多样的处理手段。多波多分量(转换波三维三分量地震勘探方法技术研究，唐建明，2010)地震勘探技术能够获取更丰富的波动信息，可以同时利用纵波、横波和转换波资料，实现地下构造异常更精确的成像结果。极化分析(王勃，2014)是研究地震波的空间质点振动问题，结合地震波入射角、传播的方位角等路径信息，在偏移成像时，主要通过分析接收点处质点振动方向与反射射线方向相关关系，完成极化成像，成像结果更加符合实际地质分布情况。将转换波技术、极化偏移技术应用到工程地震勘探资料处理中，能够实现高精度高分辨率偏移成像，有助于对现场数据的解释及工程地震勘探的发展。

技术实现要素：

本发明所要解决的技术问题在于提供了一种能够从多个角度获取更高精度和分辨率的城市工程地震探测剖面的基于多分量观测系统的城市工程地震探测综合成像方法。

本发明是通过以下技术方案解决上述技术问题的：一种基于多分量观测系统的城市工程地震探测综合成像方法，包括如下步骤：

(1)在探测区域线性等间隔布置若干接收点，每个接收点布置一个三分量检波器，包括两个水平分量和一个垂直分量，激发点布置在排列中心位置，或者在排列中间等间隔布置多个激发点，以排列的中心位置作为测线的测点位置，排列内所有激发点的探测数据都作为当前测点的采集数据，排列内的所有的激发点激发采集完成后，则当前测点数据采集结束，然后将整个排列向前移动，进行下一测点的数据采集，依次沿探测测线进行所有测点的数据采集，完成最终数据采集工作；

(2)对采集数据进行常规的数据处理，并进行速度分析，完成探测区域的速度建模；

(3)对探测区域进行偏移成像计算，采用以下计算方法中的任一种或者两者或者两者以上进行计算：单点多次覆盖偏移成像方法、反射偏移成像方法、转换波散射偏移成像方法、极化散射偏移成像方法，得出偏移成像剖面结果，基于多个剖面结果进行综合评价。

作为优选的方案，所述接收点的个数为4～16个，接收点间距控制在0.5～1米左右。

作为优选的方案，激发点采用锤击或者其他无损或微损震源。

作为优选的方案，所述步骤(1)中整个排列向前移动的距离为相邻接收点之间间距一半的倍数，移动的最大距离不大于距离最远的两个接收点之间的间距的一半。

作为优选的方案，探测区域内布置多条测线，各条测线并排平行排列。

作为优选的方案，所述步骤(3)中，偏移成像计算方法包括：

首先将排列下方探测区域进行网格剖分，沿测线方向的网格大小定义为Δx，深度方向的网格大小定义为Δz，即在X-Z平面形成二维网格点，把每一个网格点看成是一个反射点，则任意网格反射点G的反射波或绕射波旅行时为：其中，v为地震波速度，且v根据步骤(2)创建的速度模型获得，当使用转换波计算时，v按照相应的转换波类型速度计算，转换波类型为纵波时，其速度即为v，转换波为横波时，按照公式v_s＝ρv(其中：v_s为转换横波速度，ρ＝0.2～0.7)的结果作为计算时的速度。j＝1,2,3,…m，m为参与叠加的所有接收点的数量，z为网格反射点的垂直深度，t_ij为计算网格反射点处的第i炮第j个接收点的绕射波旅行时，为第i炮的坐标，为第j个接收点的坐标，x_g为网格反射点沿测线方向的坐标，把所有参与的叠加记录道上距离激发时刻t_ij的时刻的振幅值a_ij叠加到网格反射点上，作为此点的总振幅值A_i，即：依次扫描X-Z平面内的所有网格反射点，计算其叠加幅度值A_i。记录道是实际检波器接收的信号记录，在记录数据上取振幅值。

作为优选的方案，计算单点多次覆盖偏移成像方法：每个排列仅对其正下方的深度点位置进行成像，只使用接收点的垂直分量数据，首先找到X-Z平面网格反射点对应的水平面上沿测线方向的测点，然后根据测点位置找到此测点对应的排列数据集C_ij，i为排列内的激发点，j为第i激发点的排列内的所有接收点，C_ij即为参与叠加计算的所有记录道，最后按照公式进行网格反射点的振幅叠加计算。

作为优选的方案，计算反射偏移成像方法：根据测点进行所有排列的计算，单个排列内的反射区域有限，排列下方只有部分区域能够反射入射波，并由接收点检波器接收，并且只使用接收点的垂直分量数据，首先根据排列内激发点位置x_S和接收点位置x_R计算X-Z平面内网格反射点的横向位置，即然后在此横向位置处依次扫描计算网格反射点的深度z_g＝n×Δz，n为深度方向的网格，根据旅行时计算公式获得t_ij，并将距离激发时刻t_ij的时刻对应的幅度a_ij叠加到网格反射点(x_g,z_g)上，依次计算当前排列内的所有记录道，则当前测点排列的反射区域计算完成，再进行下一个测点排列的计算，直到所有测点排列全部计算，在X-Z平面得到所有网格反射点的叠加振幅值A_i，形成反射偏移成像结果。

作为优选的方案，计算转换波散射偏移成像方法：根据三分量检波器的布置条件，沿测线方向X的检波器接收方向为：L_RX:(1，0，0)，沿Y的检波器接收方向为L_RY:(0，1，0)，沿Z指向地下的检波器接收方向为L_RZ:(0，0，1)，纵波入射纵波反射的转换纵波方向为L_PP:(α，β，γ)，其中α、β、γ分别是转换纵波与坐标轴X、Y、Z的空间夹角的余弦，它们的值通过X-Z平面内已知的网格反射点位置坐标(x_g,y_g,z_g)和接收点位置坐标(x_R,y_R,z_R)进行向量计算得到，纵波入射横波SH波反射的转换SH波方向为：L_PSH:(β，-α，0)，纵波入射横波SV波反射的转换SV波方向为：L_PSV:(γα，γβ，α²-γβ)，当进行纵波入射纵波反射的转换波偏移成像计算时，绕射扫描偏移叠加的总振幅值为：其中k_pz＝L_PP·L_RZ，‘·’表示向量点乘，k_px＝L_PP·L_RX，k_py＝L_PP·L_RY。a_xij、a_yij、a_zij分别为三分量检波器接收的t_ij时刻的振幅值，当进行纵波入射横波SH反射的转换波偏移成像计算时，绕射扫描偏移叠加的总振幅值为：其中k_shz＝L_PSH·L_RZ，k_shx＝L_PSH·L_RX，k_shy＝L_PSH·L_RY，当进行纵波入射横波SV反射的转换波偏移成像计算时，绕射扫描偏移叠加的总振幅值为：其中k_shz＝L_PSV·L_RZ，k_svx＝L_PSV·L_RX，k_svy＝L_PSV·L_RY，经过上述计算，得到纵波入射纵波反射、横波SH反射和横波SV反射的转换波偏移成像结果。

作为优选的方案，计算极化散射偏移成像方法：首先计算各接收点处数据采样点的三分量数据主极化方向方位角、倾角特征参数，任意接收点处数据采样点的主极化方向为：其中θ和分别是主极化方向的方位角和倾角特征参数，主极化方向分别与X水平分量检波器、Y水平分量检波器和Z垂向分量检波器的权值计算方法为：k_mx＝L_M·L_RX，当进行纵波入射纵波反射偏移计算时，绕射扫描偏移总叠加幅度为当进行纵波入射横波SH反射偏移计算时，绕射扫描偏移总叠加幅度为当进行纵波入射横波SV反射偏移计算时，绕射扫描偏移总叠加幅度为

本发明相比现有技术具有以下优点：其利用多分量探测数据，简化转换波数据处理及极化偏移成像方法，并将其应用到成像计算当中，在实际探测区域内可以同时获取单点多次覆盖偏移成像、反射偏移成像、以及基于转换波和极化偏移技术的散射偏移成像结果，采用多种成像方法对探测区域进行综合成像，获取高精度高分辨率的城市工程地震探测剖面，便于对城市地下异常体的解释与对比，从而能反映真实的地下结构及异常位置。

附图说明

图1a为数据采集观测系统及其移动示意图。

图1b为计算空间坐标系示意图。

图2为综合成像区域示意图。

图3为单点多次覆盖偏移成像剖面结果。

图4为反射偏移成像剖面结果。

图5为转换波散射偏移成像剖面结果。

图6为极化散射偏移成像剖面结果。

具体实施方式

下面对本发明的实施例作详细说明，本实施例在以本发明技术方案为前提下进行实施，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述的实施例。

同时参考图1、图2所示，一种基于多分量观测系统的城市工程地震探测综合成像方法，其基本实施步骤包括：

(1)在探测区域线性等间隔布置若干个接收点，优选的，接收点的数量以4～16个为宜，接收点间距控制在0.5～1米左右为宜，每个接收点布置一个三分量检波器，包括两个水平分量和一个垂直分量，激发点布置在排列中心位置，采用锤击震源，也可以在排列中间等间隔布置多个激发点，以排列的中心位置作为测线的测点位置，排列内所有激发点的探测数据都作为当前测点的采集数据，排列内的所有的激发点激发采集完成后，则当前测点数据采集结束，然后将整个排列向前移动接收点间距大小一半的倍数的距离，进行下一测点的数据采集，移动的最大距离不大于距离最远的两个接收点之间的间距的一半，依次沿探测测线进行所有测点的数据采集，完成最终数据采集工作。

(2)对采集数据进行数据修饰、数字滤波等常规预处理，并进行速度分析，完成探测区域的速度建模，速度模型可以是简化的匀速度模型或者是根据共中心点道集(CMP)或共散射点道集(CSP)计算的扫描速度谱模型。例如对CMP道集来说，炮检对的旅行时间t_x是炮检距x的函数，当炮检距不很大时，近似为x的双曲线函数：其中v称为叠加速度，t₀为零偏移距旅行时，以t₀参变量为第一层循环，以v参变量为第二层循环，在第二层循环中计算CMP道集所有炮检对相应的旅行时t_x，同时取t_x对应的幅度值进行叠加作为某一v时的叠加值，扫描所有的v第二层循环结束，然后扫描所有的t₀第一层循环结束，得到速度谱图，然后拾取最大叠加值对应的v和t₀作为某CMP道集的速度模型曲线，最后对测线的多个CMP速度模型曲线进行线性插值，获得速度模型。

(3)对探测区域进行偏移成像计算。参考图2所示，观测系统排列具有一定的长度，计算不同的偏移成像时，其计算区域范围是不同的：对于单点多次覆盖偏移成像，其只计算排列中心位置对应的下方(即灰色圆圈位置)深度的叠加结果，地震波射线路径如图2中的粗黑色点线所示；对于反射偏移成像，排列内计算区域如图2中两个带有灰色填充的矩形框所示，其地震波射线路径如图2中的细点划线所示；对于散射偏移成像，其计算区域覆盖整个排列(图2中不带边框的长灰色矩形块)，地震波射线路径如图2中的细实线所示。首先将排列下方探测区域进行网格剖分，沿测线方向的网格大小定义为Δx，深度方向的网格大小定义为Δz，即在X-Z平面形成二维网格点，把每一个网格点看成是一个反射点，则任意网格反射点G的反射波或绕射波旅行时为：其中，v为地震波速度，且v根据上述步骤(2)创建的速度模型获得，当使用转换波计算时，v按照相应的转换波类型速度计算，转换波类型为纵波时，其速度即为v，转换波为横波时，按照公式v_s＝ρv(其中：v_s为转换横波速度，ρ＝0.2～0.7)的结果作为计算时的速度。j＝1,2,3,…m，m为参与叠加的所有接收点的数量，z为网格反射点的垂直深度，t_ij为计算网格反射点处的第i炮第j个接收点的绕射波旅行时，为第i炮的坐标，为第j个接收点的坐标，x_g为网格反射点沿测线方向的坐标。把所有参与叠加的记录道上距离激发时刻t_ij的时刻的振幅值a_ij叠加到网格反射点上，作为此点的总振幅值A_i，即：依次扫描X-Z平面内的所有网格反射点，计算其叠加幅度值A_i，即可得到以下四种不同叠加幅度的偏移成像剖面。

(3.1)计算单点多次覆盖偏移成像方法：每个排列仅对其正下方的深度点位置进行成像，只使用接收点的垂直分量数据，首先找到X-Z平面网格反射点对应的水平面上沿测线方向的测点，然后根据测点位置找到此测点对应的排列数据集C_ij，i为排列内的激发点，j为第i激发点的排列内的所有接收点，C_ij即为参与叠加计算的所有记录道，最后按照公式进行网格反射点的振幅叠加计算。

(3.2)计算反射偏移成像方法：根据测点进行所有排列的计算，单个排列内的反射区域有限，排列下方只有部分区域能够反射入射波，并由接收点检波器接收，并且只使用接收点的垂直分量数据，首先根据排列内激发点位置x_S和接收点位置x_R计算X-Z平面内网格反射点的横向位置，即然后在此横向位置处依次扫描计算网格反射点的深度z_g＝n×Δz，n为深度方向的网格，根据旅行时计算公式获得t_ij，并将距离激发时刻t_ij的时刻对应的幅度a_ij叠加到网格反射点(x_g,z_g)上，依次计算当前排列内的所有记录道，则当前测点排列的反射区域计算完成，再进行下一个测点排列的计算，直到所有测点排列全部计算，在X-Z平面得到所有网格反射点的叠加振幅值A_i，形成反射偏移成像结果。

(3.3)计算转换波散射偏移成像方法：根据三分量检波器的布置条件，沿测线方向X的检波器接收方向为：L_RX:(1，0，0)，沿Y的检波器接收方向为L_RY:(0，1，0)，沿Z指向地下的检波器接收方向为L_RZ:(0，0，1)。纵波入射纵波反射的转换纵波方向为L_PP:(α，β，γ)，其中α、β、γ分别是转换纵波与坐标轴X、Y、Z的空间夹角的余弦，它们的值通过X-Z平面内已知的网格反射点位置坐标(x_g,y_g,z_g)和接收点位置坐标(x_R,y_R,z_R)进行向量计算得到，纵波入射横波SH波反射的转换SH波方向为：L_PSH:(β，-α，0)，纵波入射横波SV波反射的转换SV波方向为：L_PSV:(γα，γβ，α²-γβ)。当进行纵波入射纵波反射的转换波偏移成像计算时，绕射扫描偏移叠加的总振幅值为：其中k_pz＝L_PP·L_RZ，‘·’表示向量点乘，k_px＝L_PP·L_RX，k_py＝L_PP·L_RY。a_xij、a_yij、a_zij分别为三分量检波器接收的t_ij时刻的振幅值。当进行纵波入射横波SH反射的转换波偏移成像计算时，绕射扫描偏移叠加的总振幅值为：其中k_shz＝L_PSH·L_RZ，k_shx＝L_PSH·L_RX，k_shy＝L_PSH·L_RY。当进行纵波入射横波SV反射的转换波偏移成像计算时，绕射扫描偏移叠加的总振幅值为：其中k_shz＝L_PSV·L_RZ，k_svx＝L_PSV·L_RX，k_svy＝L_PSV·L_RY，经过上述计算，可以得到纵波入射纵波反射、横波SH反射和横波SV反射的转换波偏移成像结果。

(3.4)计算极化散射偏移成像方法：首先计算各接收点处数据采样点的三分量数据主极化方向方位角、倾角特征参数。其计算步骤为：假设三分量接收检波器的记录数据分别为两个水平分量P_X(t)、P_Y(t)及垂直分量P_Z(t)，首先进行希尔伯特变换得到复地震道，HP_X(t)＝P_X(t)+j∪(P_X(t))，HP_Y(t)＝P_Y(t)+j∪(P_Y(t))，HP_Z(t)＝P_Z(t)+j∪(P_Z(t))，其中符号∪表示希尔伯特变换，j为虚数单位；然后构造复协方差矩阵，即Q(t)＝M^*(t)·M(t)，其中M(t)＝(HP_X(t),HP_Y(t),HP_Z(t))，符号*表示复共轭转置；对于协方差矩阵Q(t)存在特征方程其中λ_i(t)为某时刻对应的三个特征值，(V_i^x(t),V_i^y(t),V_i^z(t))为λ_i(t)对应的特征向量，I为单位矩阵，由于协方差矩阵Q(t)为厄米特共轭矩阵，所以λ_i(t)均为非负实数，其中最大特征值λ(t)对应的特征向量V(t)即为主极化方向向量；最后计算特征参数，对主极化向量V(t)进行归一化得到(V₁(t),V₂(t),V₃(t))，则任意时刻t对应的方位角倾角符号Re表示取复数的实部运算。特征参数计算完成后，则任意接收点处数据采样点的主极化方向为：其中θ和分别是某时刻主极化方向的方位角和倾角特征参数，主极化方向分别与X水平分量检波器、Y水平分量检波器和Z垂向分量检波器的权值计算方法为：k_mx＝L_M·L_RX，当进行纵波入射纵波反射偏移计算时，绕射扫描偏移总叠加幅度为当进行纵波入射横波SH反射偏移计算时，绕射扫描偏移总叠加幅度为当进行纵波入射横波SV反射偏移计算时，绕射扫描偏移总叠加幅度为

(4)按照以上步骤(3)的4种计算方法中的任一种或者两者或者两者以上进行计算，得出偏移成像剖面结果，对探测区域进行后期解释，多个剖面结果便于综合评价，从而能反映真实的地下结构及异常位置。

实施例一

以某地铁一段实际探测项目为例，其具体操作方式包括：

(1)三分量线性观测系统排列总共布置4个接收点，每个接收点布置一个三分量数字检波器，接收点间距0.5米，激发点在排列中心位置，在排列中心的激发点位置激发一次，完成当前测点数据的采集，然后将整个排列沿测线向前移动0.25米，进行下一个测点的数据采集，依次沿测线移动总共采集180个测点的数据采集工作。

(2)对采集的数据进行振幅恢复、数字滤波、提取道集，进而创建速度模型。

(3)分别计算探测区域内的单点多次覆盖偏移成像、反射偏移成像、转换波散射偏移成像以及极化散射偏移成像。

图3-图6为同一地点的4种成像方法的结果，可以看出图3、图4对于地下10米以浅的分辨率高，图5、图6对于地下10米以深的信噪比高；4张图可以通过地质资料对比后综合确定最终偏移成像结果。

实施例二

以某地铁一段实际探测项目为例，其具体操作方式包括：

(1)三分量线性观测系统排列总共布置8个接收点，每个接收点布置一个三分量数字检波器，接收点间距1米，排列中间等间隔布置3个激发点，在排列中间的各激发点位置各激发一次，完成当前测点数据的采集，然后将整个排列沿测线向前移动1米，进行下一个测点的数据采集，依次沿测线移动并完成60个测点的数据采集工作。

(2)对采集的数据进行振幅恢复、数字滤波、提取道集，进而创建速度模型。

(3)分别计算探测区域内的单点多次覆盖偏移成像、反射偏移成像、转换波散射偏移成像以及极化散射偏移成像。

实施例三

对要求更高精度的区域探测时，可以在探测区域内布置多条测线，各条测线对应的X轴坐标相同，Y轴坐标不同，同时在排列内等间隔布置多个激发点，包含多条测线时，可以进行三维网格划分，获取三维偏移数据体。该实施例与实施例一和实施例二的不同之处在于，实施例三增加了测线密度，对成像结果具有更高的精度。在进行三维绕射扫描偏移叠加时，则需要对探测区域进行网格体的划分，其他的计算方法和二维绕射扫描偏移叠加一致。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。