CN106443765B - 基于多分量观测系统的城市工程地震探测综合成像方法 - Google Patents
基于多分量观测系统的城市工程地震探测综合成像方法 Download PDFInfo
- Publication number
- CN106443765B CN106443765B CN201610779654.8A CN201610779654A CN106443765B CN 106443765 B CN106443765 B CN 106443765B CN 201610779654 A CN201610779654 A CN 201610779654A CN 106443765 B CN106443765 B CN 106443765B
- Authority
- CN
- China
- Prior art keywords
- wave
- point
- arrangement
- grid
- imaging method
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000003384 imaging method Methods 0.000 title claims abstract description 94
- 238000013508 migration Methods 0.000 claims abstract description 61
- 230000005012 migration Effects 0.000 claims abstract description 61
- 230000010287 polarization Effects 0.000 claims abstract description 27
- 238000012545 processing Methods 0.000 claims abstract description 8
- 238000004364 calculation method Methods 0.000 claims description 14
- 238000000034 method Methods 0.000 claims description 14
- 238000001514 detection method Methods 0.000 claims description 13
- 238000006243 chemical reaction Methods 0.000 claims description 12
- 230000005284 excitation Effects 0.000 claims description 10
- 238000000205 computational method Methods 0.000 claims description 7
- 238000005070 sampling Methods 0.000 claims description 6
- 238000004513 sizing Methods 0.000 claims description 6
- 238000004458 analytical method Methods 0.000 claims description 5
- 241001637516 Polygonia c-album Species 0.000 claims description 3
- 238000007689 inspection Methods 0.000 claims description 3
- 241000208340 Araliaceae Species 0.000 claims 1
- 235000005035 Panax pseudoginseng ssp. pseudoginseng Nutrition 0.000 claims 1
- 235000003140 Panax quinquefolius Nutrition 0.000 claims 1
- 235000008434 ginseng Nutrition 0.000 claims 1
- 238000005516 engineering process Methods 0.000 abstract description 9
- 230000002159 abnormal effect Effects 0.000 abstract description 5
- 230000002547 anomalous effect Effects 0.000 abstract description 2
- 230000002146 bilateral effect Effects 0.000 abstract 1
- 238000010276 construction Methods 0.000 description 6
- 239000011159 matrix material Substances 0.000 description 5
- 238000011161 development Methods 0.000 description 4
- 238000001914 filtration Methods 0.000 description 4
- 238000013316 zoning Methods 0.000 description 4
- 230000008569 process Effects 0.000 description 3
- 238000013480 data collection Methods 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 238000000605 extraction Methods 0.000 description 2
- 239000002245 particle Substances 0.000 description 2
- 230000000149 penetrating effect Effects 0.000 description 2
- 238000011084 recovery Methods 0.000 description 2
- 238000009412 basement excavation Methods 0.000 description 1
- 239000003245 coal Substances 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000003708 edge detection Methods 0.000 description 1
- 235000013399 edible fruits Nutrition 0.000 description 1
- 239000000686 essence Substances 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000009434 installation Methods 0.000 description 1
- 238000011835 investigation Methods 0.000 description 1
- 239000000463 material Substances 0.000 description 1
- 239000002184 metal Substances 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 239000003607 modifier Substances 0.000 description 1
- 238000013439 planning Methods 0.000 description 1
- 238000003969 polarography Methods 0.000 description 1
- 239000002689 soil Substances 0.000 description 1
- 238000001228 spectrum Methods 0.000 description 1
- 238000012876 topography Methods 0.000 description 1
- 230000017105 transposition Effects 0.000 description 1
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/02—Generating seismic energy
- G01V1/143—Generating seismic energy using mechanical driving means, e.g. motor driven shaft
- G01V1/147—Generating seismic energy using mechanical driving means, e.g. motor driven shaft using impact of dropping masses
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/16—Receiving elements for seismic signals; Arrangements or adaptations of receiving elements
- G01V1/20—Arrangements of receiving elements, e.g. geophone pattern
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/30—Analysis
- G01V1/306—Analysis for determining physical properties of the subsurface, e.g. impedance, porosity or attenuation profiles
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/32—Transforming one recording into another or one representation into another
- G01V1/325—Transforming one representation into another
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/40—Transforming data representation
- G01V2210/48—Other transforms
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/60—Analysis
- G01V2210/62—Physical property of subsurface
Landscapes
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Acoustics & Sound (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
本发明公开了一种基于多分量观测系统的城市工程地震探测综合成像方法，基本步骤包括：首先使用中间激发双边接收的观测系统完成多分量数据的采集工作；然后对采集数据进行常规数据处理并建立速度模型；最后基于创建的速度模型和基于绕射扫描偏移、转换波和极化偏移技术对探测区域进行多角度成像，获取多种分辨率的成像剖面。本发明的优点在于：采用多种成像方法对探测区域进行综合成像，获取高精度高分辨率的城市工程地震探测剖面，便于对城市地下异常体的解释与对比，从而能反映真实的地下结构及异常位置。
Description
技术领域
本发明属于工程地震数据成像处理领域，具体涉及一种城市工程地震探测的成像方法。
背景技术
随着我国国民经济建设的飞速发展，城市化进程不断深入，城市工程建设在规划、设计、施工阶段都必须对建设区域内的地质情况及地下埋设物有一个系统的勘察，在建设工程中及建成后还必须对工程质量进行检测和监测。工程物探是用于探查工程地质问题(如岩土分层、基岩面形态、岩溶分布、断层裂隙、孤石及复杂地基土等)、环境地质问题(如地下采空塌陷、地面沉降、废弃箱涵、暗河与管线等)的有效手段，其在城市工程建设的作用越来越重要，特别是工程地震勘探技术，由于其采用弹性波方法，具有抗金属与电器干扰能力强，探测深度较大的显著特征，是目前城市工程物探中急需发展的方法。
工程地震勘探在城市地下探测的应用，对工程地震数据处理的精度提出了更高的要求，其数据处理结果将直接影响后期资料解释的准确性和可靠度。目前，工程地震勘探的方法主要是建立在野外石油、煤炭等能源地震勘探的基础上，其数据处理一般包括前期的预处理、数字滤波、速度分析、叠加偏移等，特别是地震偏移技术，它是用计算机按一定的计算方法对观测数据进行处理，使之成为反映地下地质分层界面位置及反射系数值的反射界面的像，是地震勘探数据处理的关键技术，其在一定程度上代表了地震勘探技术的发展程度。绕射扫描偏移是建立在射线偏移基础上使反射波自动归位到真实位置上的一种方法，根据惠更斯原理，地下每一个反射点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，转换波为横波时，按照公式vs＝ρv(其中：vs为转换横波速度，ρ＝0.2～0.7)的结果作为计算时的速度。j＝1,2,3,…m，m为参与叠加的所有接收点的数量，z为网格反射点的垂直深度，tij为计算网格反射点处的第i炮第j个接收点的绕射波旅行时，为第i炮的坐标，为第j个接收点的坐标，xg为网格反射点沿测线方向的坐标，把所有参与的叠加记录道上距离激发时刻tij的时刻的振幅值aij叠加到网格反射点上，作为此点的总振幅值Ai，即：依次扫描X-Z平面内的所有网格反射点，计算其叠加幅度值Ai。记录道是实际检波器接收的信号记录，在记录数据上取振幅值。
作为优选的方案，计算单点多次覆盖偏移成像方法：每个排列仅对其正下方的深度点位置进行成像，只使用接收点的垂直分量数据，首先找到X-Z平面网格反射点对应的水平面上沿测线方向的测点，然后根据测点位置找到此测点对应的排列数据集Cij，i为排列内的激发点，j为第i激发点的排列内的所有接收点，Cij即为参与叠加计算的所有记录道，最后按照公式进行网格反射点的振幅叠加计算。
作为优选的方案，计算反射偏移成像方法：根据测点进行所有排列的计算，单个排列内的反射区域有限，排列下方只有部分区域能够反射入射波，并由接收点检波器接收，并且只使用接收点的垂直分量数据，首先根据排列内激发点位置xS和接收点位置xR计算X-Z平面内网格反射点的横向位置，即然后在此横向位置处依次扫描计算网格反射点的深度zg＝n×Δz，n为深度方向的网格，根据旅行时计算公式获得tij，并将距离激发时刻tij的时刻对应的幅度aij叠加到网格反射点(xg,zg)上，依次计算当前排列内的所有记录道，则当前测点排列的反射区域计算完成，再进行下一个测点排列的计算，直到所有测点排列全部计算，在X-Z平面得到所有网格反射点的叠加振幅值Ai，形成反射偏移成像结果。
作为优选的方案，计算转换波散射偏移成像方法：根据三分量检波器的布置条件，沿测线方向X的检波器接收方向为：LRX:(1，0，0)，沿Y的检波器接收方向为LRY:(0，1，0)，沿Z指向地下的检波器接收方向为LRZ:(0，0，1)，纵波入射纵波反射的转换纵波方向为LPP:(α，β，γ)，其中α、β、γ分别是转换纵波与坐标轴X、Y、Z的空间夹角的余弦，它们的值通过X-Z平面内已知的网格反射点位置坐标(xg,yg,zg)和接收点位置坐标(xR,yR,zR)进行向量计算得到，纵波入射横波SH波反射的转换SH波方向为：LPSH:(β，-α，0)，纵波入射横波SV波反射的转换SV波方向为：LPSV:(γα，γβ，α2-γβ)，当进行纵波入射纵波反射的转换波偏移成像计算时，绕射扫描偏移叠加的总振幅值为：其中kpz＝LPP·LRZ，‘·’表示向量点乘，kpx＝LPP·LRX，kpy＝LPP·LRY。axij、ayij、azij分别为三分量检波器接收的tij时刻的振幅值，当进行纵波入射横波SH反射的转换波偏移成像计算时，绕射扫描偏移叠加的总振幅值为：其中kshz＝LPSH·LRZ，kshx＝LPSH·LRX，kshy＝LPSH·LRY，当进行纵波入射横波SV反射的转换波偏移成像计算时，绕射扫描偏移叠加的总振幅值为：其中kshz＝LPSV·LRZ，ksvx＝LPSV·LRX，ksvy＝LPSV·LRY，经过上述计算，得到纵波入射纵波反射、横波SH反射和横波SV反射的转换波偏移成像结果。
作为优选的方案，计算极化散射偏移成像方法：首先计算各接收点处数据采样点的三分量数据主极化方向方位角、倾角特征参数，任意接收点处数据采样点的主极化方向为：其中θ和分别是主极化方向的方位角和倾角特征参数，主极化方向分别与X水平分量检波器、Y水平分量检波器和Z垂向分量检波器的权值计算方法为：kmx＝LM·LRX，当进行纵波入射纵波反射偏移计算时，绕射扫描偏移总叠加幅度为当进行纵波入射横波SH反射偏移计算时，绕射扫描偏移总叠加幅度为当进行纵波入射横波SV反射偏移计算时，绕射扫描偏移总叠加幅度为
本发明相比现有技术具有以下优点：其利用多分量探测数据，简化转换波数据处理及极化偏移成像方法，并将其应用到成像计算当中，在实际探测区域内可以同时获取单点多次覆盖偏移成像、反射偏移成像、以及基于转换波和极化偏移技术的散射偏移成像结果，采用多种成像方法对探测区域进行综合成像，获取高精度高分辨率的城市工程地震探测剖面，便于对城市地下异常体的解释与对比，从而能反映真实的地下结构及异常位置。
附图说明
图1a为数据采集观测系统及其移动示意图。
图1b为计算空间坐标系示意图。
图2为综合成像区域示意图。
图3为单点多次覆盖偏移成像剖面结果。
图4为反射偏移成像剖面结果。
图5为转换波散射偏移成像剖面结果。
图6为极化散射偏移成像剖面结果。
具体实施方式
下面对本发明的实施例作详细说明，本实施例在以本发明技术方案为前提下进行实施，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述的实施例。
同时参考图1、图2所示，一种基于多分量观测系统的城市工程地震探测综合成像方法，其基本实施步骤包括：
(1)在探测区域线性等间隔布置若干个接收点，优选的，接收点的数量以4～16个为宜，接收点间距控制在0.5～1米左右为宜，每个接收点布置一个三分量检波器，包括两个水平分量和一个垂直分量，激发点布置在排列中心位置，采用锤击震源，也可以在排列中间等间隔布置多个激发点，以排列的中心位置作为测线的测点位置，排列内所有激发点的探测数据都作为当前测点的采集数据，排列内的所有的激发点激发采集完成后，则当前测点数据采集结束，然后将整个排列向前移动接收点间距大小一半的倍数的距离，进行下一测点的数据采集，移动的最大距离不大于距离最远的两个接收点之间的间距的一半，依次沿探测测线进行所有测点的数据采集，完成最终数据采集工作。
(2)对采集数据进行数据修饰、数字滤波等常规预处理，并进行速度分析，完成探测区域的速度建模，速度模型可以是简化的匀速度模型或者是根据共中心点道集(CMP)或共散射点道集(CSP)计算的扫描速度谱模型。例如对CMP道集来说，炮检对的旅行时间tx是炮检距x的函数，当炮检距不很大时，近似为x的双曲线函数：其中v称为叠加速度，t0为零偏移距旅行时，以t0参变量为第一层循环，以v参变量为第二层循环，在第二层循环中计算CMP道集所有炮检对相应的旅行时tx，同时取tx对应的幅度值进行叠加作为某一v时的叠加值，扫描所有的v第二层循环结束，然后扫描所有的t0第一层循环结束，得到速度谱图，然后拾取最大叠加值对应的v和t0作为某CMP道集的速度模型曲线，最后对测线的多个CMP速度模型曲线进行线性插值，获得速度模型。
(3)对探测区域进行偏移成像计算。参考图2所示，观测系统排列具有一定的长度，计算不同的偏移成像时，其计算区域范围是不同的：对于单点多次覆盖偏移成像，其只计算排列中心位置对应的下方(即灰色圆圈位置)深度的叠加结果，地震波射线路径如图2中的粗黑色点线所示；对于反射偏移成像，排列内计算区域如图2中两个带有灰色填充的矩形框所示，其地震波射线路径如图2中的细点划线所示；对于散射偏移成像，其计算区域覆盖整个排列(图2中不带边框的长灰色矩形块)，地震波射线路径如图2中的细实线所示。首先将排列下方探测区域进行网格剖分，沿测线方向的网格大小定义为Δx，深度方向的网格大小定义为Δz，即在X-Z平面形成二维网格点，把每一个网格点看成是一个反射点，则任意网格反射点G的反射波或绕射波旅行时为：其中，v为地震波速度，且v根据上述步骤(2)创建的速度模型获得，当使用转换波计算时，v按照相应的转换波类型速度计算，转换波类型为纵波时，其速度即为v，转换波为横波时，按照公式vs＝ρv(其中：vs为转换横波速度，ρ＝0.2～0.7)的结果作为计算时的速度。j＝1,2,3,…m，m为参与叠加的所有接收点的数量，z为网格反射点的垂直深度，tij为计算网格反射点处的第i炮第j个接收点的绕射波旅行时，为第i炮的坐标，为第j个接收点的坐标，xg为网格反射点沿测线方向的坐标。把所有参与叠加的记录道上距离激发时刻tij的时刻的振幅值aij叠加到网格反射点上，作为此点的总振幅值Ai，即：依次扫描X-Z平面内的所有网格反射点，计算其叠加幅度值Ai，即可得到以下四种不同叠加幅度的偏移成像剖面。
(3.1)计算单点多次覆盖偏移成像方法：每个排列仅对其正下方的深度点位置进行成像，只使用接收点的垂直分量数据，首先找到X-Z平面网格反射点对应的水平面上沿测线方向的测点，然后根据测点位置找到此测点对应的排列数据集Cij，i为排列内的激发点，j为第i激发点的排列内的所有接收点，Cij即为参与叠加计算的所有记录道，最后按照公式进行网格反射点的振幅叠加计算。
(3.2)计算反射偏移成像方法：根据测点进行所有排列的计算，单个排列内的反射区域有限，排列下方只有部分区域能够反射入射波，并由接收点检波器接收，并且只使用接收点的垂直分量数据，首先根据排列内激发点位置xS和接收点位置xR计算X-Z平面内网格反射点的横向位置，即然后在此横向位置处依次扫描计算网格反射点的深度zg＝n×Δz，n为深度方向的网格，根据旅行时计算公式获得tij，并将距离激发时刻tij的时刻对应的幅度aij叠加到网格反射点(xg,zg)上，依次计算当前排列内的所有记录道，则当前测点排列的反射区域计算完成，再进行下一个测点排列的计算，直到所有测点排列全部计算，在X-Z平面得到所有网格反射点的叠加振幅值Ai，形成反射偏移成像结果。
(3.3)计算转换波散射偏移成像方法：根据三分量检波器的布置条件，沿测线方向X的检波器接收方向为：LRX:(1，0，0)，沿Y的检波器接收方向为LRY:(0，1，0)，沿Z指向地下的检波器接收方向为LRZ:(0，0，1)。纵波入射纵波反射的转换纵波方向为LPP:(α，β，γ)，其中α、β、γ分别是转换纵波与坐标轴X、Y、Z的空间夹角的余弦，它们的值通过X-Z平面内已知的网格反射点位置坐标(xg,yg,zg)和接收点位置坐标(xR,yR,zR)进行向量计算得到，纵波入射横波SH波反射的转换SH波方向为：LPSH:(β，-α，0)，纵波入射横波SV波反射的转换SV波方向为：LPSV:(γα，γβ，α2-γβ)。当进行纵波入射纵波反射的转换波偏移成像计算时，绕射扫描偏移叠加的总振幅值为：其中kpz＝LPP·LRZ，‘·’表示向量点乘，kpx＝LPP·LRX，kpy＝LPP·LRY。axij、ayij、azij分别为三分量检波器接收的tij时刻的振幅值。当进行纵波入射横波SH反射的转换波偏移成像计算时，绕射扫描偏移叠加的总振幅值为：其中kshz＝LPSH·LRZ，kshx＝LPSH·LRX，kshy＝LPSH·LRY。当进行纵波入射横波SV反射的转换波偏移成像计算时，绕射扫描偏移叠加的总振幅值为：其中kshz＝LPSV·LRZ，ksvx＝LPSV·LRX，ksvy＝LPSV·LRY，经过上述计算，可以得到纵波入射纵波反射、横波SH反射和横波SV反射的转换波偏移成像结果。
(3.4)计算极化散射偏移成像方法：首先计算各接收点处数据采样点的三分量数据主极化方向方位角、倾角特征参数。其计算步骤为：假设三分量接收检波器的记录数据分别为两个水平分量PX(t)、PY(t)及垂直分量PZ(t)，首先进行希尔伯特变换得到复地震道，HPX(t)＝PX(t)+j∪(PX(t))，HPY(t)＝PY(t)+j∪(PY(t))，HPZ(t)＝PZ(t)+j∪(PZ(t))，其中符号∪表示希尔伯特变换，j为虚数单位；然后构造复协方差矩阵，即Q(t)＝M*(t)·M(t)，其中M(t)＝(HPX(t),HPY(t),HPZ(t))，符号*表示复共轭转置；对于协方差矩阵Q(t)存在特征方程其中λi(t)为某时刻对应的三个特征值，(Vi x(t),Vi y(t),Vi z(t))为λi(t)对应的特征向量，I为单位矩阵，由于协方差矩阵Q(t)为厄米特共轭矩阵，所以λi(t)均为非负实数，其中最大特征值λ(t)对应的特征向量V(t)即为主极化方向向量；最后计算特征参数，对主极化向量V(t)进行归一化得到(V1(t),V2(t),V3(t))，则任意时刻t对应的方位角倾角符号Re表示取复数的实部运算。特征参数计算完成后，则任意接收点处数据采样点的主极化方向为：其中θ和分别是某时刻主极化方向的方位角和倾角特征参数，主极化方向分别与X水平分量检波器、Y水平分量检波器和Z垂向分量检波器的权值计算方法为：kmx＝LM·LRX， 当进行纵波入射纵波反射偏移计算时，绕射扫描偏移总叠加幅度为当进行纵波入射横波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轴坐标不同，同时在排列内等间隔布置多个激发点，包含多条测线时，可以进行三维网格划分，获取三维偏移数据体。该实施例与实施例一和实施例二的不同之处在于，实施例三增加了测线密度，对成像结果具有更高的精度。在进行三维绕射扫描偏移叠加时，则需要对探测区域进行网格体的划分，其他的计算方法和二维绕射扫描偏移叠加一致。
以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。
Claims (9)
1.一种基于多分量观测系统的城市工程地震探测综合成像方法，其特征在于，包括如下步骤：
(1)在探测区域线性等间隔布置若干接收点，每个接收点布置一个三分量检波器，包括两个水平分量和一个垂直分量，激发点布置在排列中心位置，或者在排列中间等间隔布置多个激发点，以排列的中心位置作为测线的测点位置，排列内所有激发点的探测数据都作为当前测点的采集数据，排列内的所有的激发点激发采集完成后，则当前测点数据采集结束，然后将整个排列向前移动，进行下一测点的数据采集，依次沿探测测线进行所有测点的数据采集，完成最终数据采集工作，所述整个排列向前移动的距离为相邻接收点之间间距一半的倍数，移动的最大距离不大于距离最远的两个接收点之间的间距的一半；
(2)对采集数据进行常规的数据处理，并进行速度分析，完成探测区域的速度建模；
(3)对探测区域进行偏移成像计算，采用以下计算方法中的任一种或者两者或者两者以上进行计算：单点多次覆盖偏移成像方法、反射偏移成像方法、转换波散射偏移成像方法、极化散射偏移成像方法，得出偏移成像剖面结果，基于多个剖面结果进行综合评价。
2.根据权利要求1所述的基于多分量观测系统的城市工程地震探测综合成像方法，其特征在于，所述接收点的个数为4～16个，接收点间距控制在0.5～1米。
3.根据权利要求1所述的基于多分量观测系统的城市工程地震探测综合成像方法，其特征在于，激发点采用锤击或者其他无损或微损震源。
4.根据权利要求1所述的基于多分量观测系统的城市工程地震探测综合成像方法，其特征在于，探测区域内布置多条测线，各条测线并排平行排列。
5.根据权利要求1所述的基于多分量观测系统的城市工程地震探测综合成像方法，其特征在于，所述步骤(3)中，偏移成像计算方法包括：
首先将排列下方探测区域进行网格剖分，沿测线方向的网格大小定义为Δx，深度方向的网格大小定义为Δz，即在X-Z平面形成二维网格点，把每一个网格点看成是一个反射点，则任意网格反射点G的反射波或绕射波旅行时为：其中，v为地震波速度，且v根据步骤(2)创建的速度模型获得，当使用转换波计算时，v按照相应的转换波类型速度计算，转换波类型为纵波时，其速度即为v，转换波为横波时，按照公式vs＝ρv的结果作为计算时的速度，其中：vs为转换横波速度，ρ＝0.2～0.7，j＝1,2,3,…m，m为参与叠加的所有接收点的数量，z为网格反射点的垂直深度，tij为计算网格反射点处的第i炮第j个接收点的绕射波旅行时，为第i炮的坐标，为第j个接收点的坐标，xg为网格反射点沿测线方向的坐标，把所有参与叠加的记录道上距离激发时刻tij的时刻的振幅值aij叠加到网格反射点上，作为此点的总振幅值Ai，即：依次扫描X-Z平面内的所有网格反射点，计算其叠加幅度值Ai。
6.根据权利要求5所述的基于多分量观测系统的城市工程地震探测综合成像方法，其特征在于，计算单点多次覆盖偏移成像方法：每个排列仅对其正下方的深度点位置进行成像，只使用接收点的垂直分量数据，首先找到X-Z平面网格反射点对应的水平面上沿测线方向的测点，然后根据测点位置找到此测点对应的排列数据集Cij，i为排列内的激发点，j为第i激发点的排列内的所有接收点，Cij即为参与叠加计算的所有记录道，最后按照公式进行网格反射点的振幅叠加计算。
7.根据权利要求5所述的基于多分量观测系统的城市工程地震探测综合成像方法，其特征在于，计算反射偏移成像方法：根据测点进行所有排列的计算，单个排列内的反射区域有限，排列下方只有部分区域能够反射入射波，并由接收点检波器接收，并且只使用接收点的垂直分量数据，首先根据排列内激发点位置xS和接收点位置xR计算X-Z平面内网格反射点的横向位置，即然后在此横向位置处依次扫描计算网格反射点的深度zg＝n×Δz，n为深度方向的网格，根据旅行时计算公式获得tij，并将距离激发时刻tij的时刻对应的幅度aij叠加到网格反射点(xg,zg)上，依次计算当前排列内的所有记录道，则当前测点排列的反射区域计算完成，再进行下一个测点排列的计算，直到所有测点排列全部计算，在X-Z平面得到所有网格反射点的叠加振幅值Ai，形成反射偏移成像结果。
8.根据权利要求5所述的基于多分量观测系统的城市工程地震探测综合成像方法，其特征在于，计算转换波散射偏移成像方法：根据三分量检波器的布置条件，沿测线方向X的检波器接收方向为：LRX:(1，0，0)，沿Y的检波器接收方向为LRY:(0，1，0)，沿Z指向地下的检波器接收方向为LRZ:(0，0，1)，纵波入射纵波反射的转换纵波方向为LPP:(α，β，γ)，其中α、β、γ分别是转换纵波与坐标轴X、Y、Z的空间夹角的余弦，它们的值通过X-Z平面内已知的网格反射点位置坐标(xg,yg,zg)和接收点位置坐标(xR,yR,zR)进行向量计算得到，纵波入射横波SH波反射的转换SH波方向为：LPSH:(β，-α，0)，纵波入射横波SV波反射的转换SV波方向为：LPSV:(γα，γβ，α2-γβ)，当进行纵波入射纵波反射的转换波偏移成像计算时，绕射扫描偏移叠加的总振幅值为：其中kpz＝LPP·LRZ，‘·’表示向量点乘，kpx＝LPP·LRX，kpy＝LPP·LRY。axij、ayij、azij分别为三分量检波器接收的tij时刻的振幅值，当进行纵波入射横波SH反射的转换波偏移成像计算时，绕射扫描偏移叠加的总振幅值为：其中kshz＝LPSH·LRZ，kshx＝LPSH·LRX，kshy＝LPSH·LRY，当进行纵波入射横波SV反射的转换波偏移成像计算时，绕射扫描偏移叠加的总振幅值为：其中kshz＝LPSV·LRZ，ksvx＝LPSV·LRX，ksvy＝LPSV·LRY，经过上述计算，得到纵波入射纵波反射、横波SH反射和横波SV反射的转换波偏移成像结果。
9.根据权利要求5所述的基于多分量观测系统的城市工程地震探测综合成像方法，其特征在于，计算极化散射偏移成像方法：首先计算各接收点处数据采样点的三分量数据主极化方向方位角、倾角特征参数，任意接收点处数据采样点的主极化方向为：其中θ和分别是主极化方向的方位角和倾角特征参数，主极化方向分别与X水平分量检波器、Y水平分量检波器和Z垂向分量检波器的权值计算方法为：kmx＝LM·LRX，当进行纵波入射纵波反射偏移计算时，绕射扫描偏移总叠加幅度为当进行纵波入射横波SH反射偏移计算时，绕射扫描偏移总叠加幅度为当进行纵波入射横波SV反射偏移计算时，绕射扫描偏移总叠加幅度为
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201610779654.8A CN106443765B (zh) | 2016-08-30 | 2016-08-30 | 基于多分量观测系统的城市工程地震探测综合成像方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201610779654.8A CN106443765B (zh) | 2016-08-30 | 2016-08-30 | 基于多分量观测系统的城市工程地震探测综合成像方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN106443765A CN106443765A (zh) | 2017-02-22 |
| CN106443765B true CN106443765B (zh) | 2018-08-28 |
Family
ID=58090361
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201610779654.8A Active CN106443765B (zh) | 2016-08-30 | 2016-08-30 | 基于多分量观测系统的城市工程地震探测综合成像方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN106443765B (zh) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111103620A (zh) * | 2019-11-20 | 2020-05-05 | 李志勇 | 一种岩巷超前探测三维偏移成像方法 |
Families Citing this family (12)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108398487A (zh) * | 2018-02-01 | 2018-08-14 | 河海大学 | 一种基于阵列弹性波技术混凝土缺陷检测方法 |
| CN108732618A (zh) * | 2018-04-23 | 2018-11-02 | 太原理工大学 | 一种用于地震监视记录上快速识别横波的装置和方法 |
| CN109375251B (zh) * | 2018-09-29 | 2021-04-13 | 山东大学 | 利用城市既有地下空间与地表的探测方法及系统 |
| CN110531417B (zh) * | 2019-08-21 | 2020-12-29 | 中国矿业大学 | 一种基于极化偏移的超前多层速度精细建模方法 |
| CN110531418B (zh) * | 2019-08-21 | 2020-11-20 | 徐州工程学院 | 一种基于希尔伯特极化成像的断点三维精细定位方法 |
| CN110780344A (zh) * | 2019-11-04 | 2020-02-11 | 北京化工大学 | 浅地表构造成像方法及装置 |
| CN111142160B (zh) * | 2019-12-30 | 2022-05-03 | 长江勘测规划设计研究有限责任公司 | 时间推移地震观测数据的分析方法及装置 |
| CN111142151B (zh) * | 2019-12-30 | 2022-05-03 | 长江勘测规划设计研究有限责任公司 | 时间推移地震观测方法及装置 |
| CN112130207B (zh) * | 2020-09-25 | 2021-07-20 | 中国科学院武汉岩土力学研究所 | 一种基于球形装药条件下由地面振动计算地下振动的方法 |
| CN112505749B (zh) * | 2020-10-19 | 2024-04-26 | 中国地质调查局南京地质调查中心(华东地质科技创新中心) | 一种基于线形台阵多次覆盖的微动数据采集方法 |
| CN112379403B (zh) * | 2020-12-14 | 2024-01-16 | 北京华晖探测科技股份有限公司 | 一种地下采空区的探测方法及系统 |
| CN114185082B (zh) * | 2021-12-02 | 2023-04-21 | 中国矿业大学 | 一种基于工作面透射地震观测的煤层下伏陷落柱探测方法 |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104570110A (zh) * | 2013-10-29 | 2015-04-29 | 中国石油化工股份有限公司 | 一种基于纵横波匹配的多分量资料联合速度分析方法 |
| CN105652320A (zh) * | 2015-12-30 | 2016-06-08 | 中国石油天然气集团公司 | 逆时偏移成像方法和装置 |
Family Cites Families (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101403797B (zh) * | 2008-11-14 | 2011-03-23 | 北京市市政工程研究院 | 地下工程施工超前地质预报系统及其预报方法 |
| CN101915938B (zh) * | 2010-07-05 | 2012-02-29 | 中国科学院地质与地球物理研究所 | 一种转换波的偏移成像方法及装置 |
| MX2014000712A (es) * | 2011-07-19 | 2014-02-20 | Halliburton Energy Serv Inc | Sistema y metodo para generacion de imagen de migracion de tensor momento. |
| CN103675908A (zh) * | 2012-09-21 | 2014-03-26 | 中国石油化工股份有限公司 | 一种海量数据图形处理器的波动方程逆时偏移成像方法 |
| CN104181581B (zh) * | 2014-08-26 | 2017-05-10 | 北京市市政工程研究院 | 基于任意排布的地震波地下工程空间观测的系统及方法 |
| CN105353406B (zh) * | 2015-10-23 | 2017-12-05 | 中国石油天然气集团公司 | 一种生成角道集的方法和装置 |
-
2016
- 2016-08-30 CN CN201610779654.8A patent/CN106443765B/zh active Active
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104570110A (zh) * | 2013-10-29 | 2015-04-29 | 中国石油化工股份有限公司 | 一种基于纵横波匹配的多分量资料联合速度分析方法 |
| CN105652320A (zh) * | 2015-12-30 | 2016-06-08 | 中国石油天然气集团公司 | 逆时偏移成像方法和装置 |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111103620A (zh) * | 2019-11-20 | 2020-05-05 | 李志勇 | 一种岩巷超前探测三维偏移成像方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN106443765A (zh) | 2017-02-22 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN106443765B (zh) | 基于多分量观测系统的城市工程地震探测综合成像方法 | |
| Forte et al. | Imaging and characterization of a carbonate hydrocarbon reservoir analogue using GPR attributes | |
| CN103645503B (zh) | 一种三维时间域照明分析及振幅补偿方法 | |
| Mordvinova et al. | Velocity structure of the lithosphere on the 2003 Mongolian-Baikal transect from SV waves | |
| Baglari et al. | A state-of-the-art review of passive MASW survey for subsurface profiling | |
| Vanlı Senkaya et al. | Integrated shallow seismic imaging of a settlement located in a historical landslide area | |
| Liu et al. | High-resolution seismic reflection survey crossing the Insubric Line into the Ivrea-Verbano Zone: Novel approaches for interpreting the seismic response of steeply dipping structures | |
| CN111103620A (zh) | 一种岩巷超前探测三维偏移成像方法 | |
| Fang et al. | Directional sensitivity of das and its effect on Rayleigh‐wave tomography: A case study in Oxnard, California | |
| CN108919351A (zh) | 基于逆时聚焦原理进行观测系统双向聚焦性的评价方法 | |
| Carrier et al. | Affordable gravity prospection calibrated on improved time-to-depth conversion of old seismic profiles for exploration of geothermal resources | |
| Fang et al. | Improving seismic remote sensing of typhoon with a three-dimensional Earth model | |
| ÇAKIR et al. | Dispersion of Rayleigh surface waves and electrical resistivities utilized to invert near surface structural heterogeneities | |
| CN103513279A (zh) | 一种基于地震波波动方程的照明分析计算方法及计算装置 | |
| Nie et al. | Integrated ERT, seismic, and electrical resistivity imaging for geological prospecting on Metro Line R3 in Qingdao, China | |
| Kolínský et al. | Velocity model of the Hronov-Poříčí Fault Zone from Rayleigh wave dispersion | |
| Han et al. | Gaussian beam summation migration of deep reflection seismic data: Numerical examples | |
| Mohamed et al. | Near-surface site characterization at Quriyat City, Sultanate of Oman using HVSR and MASW techniques | |
| Lehmann et al. | Exploration of tunnel alignment using geophysical methods to increase safety for planning and minimizing risk | |
| Bianchi et al. | A new seismic data set on the depth of the Moho in the Alps | |
| Cattaneo et al. | Propagation anomalies in Northwestern Italy by inversion of teleseismic residuals | |
| Reddy | Historical development of seismic imaging technique–an overview | |
| Botter et al. | Seismic attribute analysis of a fault zone in the Thebe field, Northwest shelf, Australia | |
| Joh et al. | Two-dimensional imaging of soil–bedrock interface by short-array beamforming technique | |
| Yue et al. | Identification and extraction of lateral target signals in tunnel geological prediction with the Karhunen-Loéve beamforming method |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| C10 | Entry into substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |