一种利用二维视电阻率数据识别地质异常体的方法

摘要

本发明涉及一种利用二维视电阻率数据识别地质异常体的方法，属于工程地球物理勘探技术领域。具体是对研究区域内用物探方法获得的视电阻率数据，利用离散二进小波变换系数的模极大值的位置同信号的局部奇异性密切相关性，辨识出信号的局部奇异点，定量分析出异常地质体的边界信息；同时，利用直方图均衡化方法增强视电阻率图像的细节对比，将二者结合进行数据分析，获取地下地质异常体的信息。该方法能使获得的图像信息更为丰富和直观，减少人为因素的影响，提高不良地质体的边界分辨精度。本发明能够对物探数据进行有效分析，获得对断层、破碎带、溶洞等异常区域更为精细化的描述，为设计、施工提供科学依据。

权利要求

1.一种利用二维视电阻率数据识别地质异常体的方法，其特征在于：基于二维视电阻率数据利用小波变换和图像增强表达相结合的数据挖掘方法进行地质异常体识别，其识别方法按以下步骤进行：①利用灰度变换函数T(r)将研究区域内用物探方法获得的任一测线上所有测点的二维视电阻率数据数组a[x][y]进行直方图均衡化，均衡化的二维视电阻率数据为s[x][y]；

其中，x为测点桩号，y为该测点下某一深度范围内任一点的深度；

②利用小波变换检测二维视电阻率数据的奇异性：

a)将任意桩号位置xp随深度变化的视电阻率数据a[xp][y]看成是随着深度y变化的一维信号f(y)；

其中，p＝1,2,3…m；

b)设f(y)为输入信号，θ(y)是某一起平滑作用的低通平滑函数，ψ(y)为θ(y)的一阶导数 记 则对尺度为2j的小波变换为：其中，j为小波变换的阶数；

c)利用Matlab软件的Findpeaks函数找出小波变换 的模极大值，确定第j阶的模极大值点的位置序列其中，设第j阶有n个极大值点，k＝1,2,3…n；

d)对第j阶的每个极大值点位置 向下在相邻尺度j-1阶上寻找所对应的极大值传播点位置 然后再以j-1阶的 为模极大值点位置再次向下寻找j-2阶上所对应的极大值传播点位置 继续向下回溯，直到寻找到第一阶上的 为止；得到的第一阶上 的位置即为第j阶的模极大值点的位置序列 所对应的异常分界点所在的位置；

e)对其他所有桩号的视电阻率数据重复b)～d)步骤执行，直到所有桩号视电阻率数据均分析完成；

③将①中直方图均衡化后的视电阻率数据图像与②中利用小波变换检测出的奇异点叠加显示，如果奇异点成团、成线状分布，并且视电阻率图像也呈明显低阻特征，可据此勾勒出异常地质体边界。

说明书

一种利用二维视电阻率数据识别地质异常体的方法

技术领域

[0001] 本发明涉及一种利用二维视电阻率数据识别地质异常体的方法，该方法用于大地电磁法物探数据的分析处理，具体是通过对大地电磁法采集的物探原始数据进行数据分析，获取地下地质异常体的信息。属于工程地球物理勘探技术领域。

背景技术

[0002] 大地电磁法是工程物探中常采用的一种物探方法，主要用于探测地下岩层中存在的不良地质信息。该方法获取的数据为二维视电阻率数据。常规物探解释基于二维视电阻率数据识别地质异常体主要根据视电阻率值的相对值、等值线形态等结合钻孔资料、地质资料等进行推断与解释，如林钧在《利用高频大地电磁法对三峡某隧道不良地质体的勘探》研究中，利用EH-4系统对视电阻率成像，结合钻孔地质资料，以及电阻率值的相对值、等值线的形态、等值线的变化梯度等综合因素来进行地质推断解释，判断地层产状、断层及岩溶发育带等。但常规物探解释受地形起伏影响较大，对解译人员经验的依赖程度较高，推断的异常地质体位置比较概略，边界难以界定，如张超在《高密度电测深法在工程勘察中的应用》中指出高密度电测深法中各岩层视电阻率受地形起伏影响较大，地表不均匀性会造成局部电性差异偏大，解释人员一定要到野外实测现场并认真做好地形地貌及地质概况记录作为解释资料时必要的参考依据。而且，根据经验分析判断往往与分析者的个体因素、经历有着很大的关联，难以做到客观、统一。

发明内容

[0003] 针对目前利用大地电磁法获取的二维视电阻率数据分析不良地质体方法存在的问题，本发明了提供了一种利用二维视电阻率数据识别地质异常体方法，其技术解决方案为：

[0004] 一种利用二维视电阻率数据识别地质异常体的方法，基于二维视电阻率数据利用小波变换和图像增强表达相结合的数据挖掘方法进行地质异常体识别，其识别方法按以下步骤进行：

[0005] ①利用灰度变换函数T(r)将研究区域内用物探方法获得的任一测线上所有测点的二维视电阻率数据数组a[x][y]进行直方图均衡化，均衡化的二维视电阻率数据为s[x][y]；

[0006] 其中，x为测点桩号，y为该测点下某一深度范围内任一点的深度。

[0007] ②利用小波变换检测二维视电阻率数据的奇异性：

[0008] a)将任意桩号位置xp随深度变化的视电阻率数据a[xp][y]看成是随着深度y变化的一维信号f(y)；

[0009] 其中，p＝1,2,3…m。

[0010] b)设f(y)为输入信号，θ(y)是某一起平滑作用的低通平滑函数，ψ(y)为θ(y)的一阶导数 记 则对尺度为2j的小波变换为：

[0011]

[0012] 其中，j为小波变换的阶数。

[0013] c)利用Matlab软件的Findpeaks函数找出小波变换 的模极大值，确定第j阶的模极大值点的位置序列

[0014] 其中，设第j阶有n个极大值点，k＝1,2,3…n。

[0015] d)对第j阶的每个极大值点 向下在相邻尺度j-1阶上寻找所对应的极大值传播点位置  然后再以j-1阶的 为模极大值点的位置序列再次向下寻找j-2阶上所对应的极大值传播点位置 继续向下回溯，直到寻找到第一阶上的 为止。得到的第一阶上的位置即为第j阶的模极大值点的位置序列 所对应的异常分界点所在的位置。

[0016] e)对其他所有桩号的视电阻率数据重复b)～d)步骤执行，直到所有桩号视电阻率数据均分析完成。

[0017] ③将①中直方图均衡化后的视电阻率数据图像与②中利用小波变换检测出的奇异点叠加显示，如果奇异点成团、成线状分布，并且视电阻率图像也呈明显低阻特征，可据此勾勒出异常地质体边界。

[0018] 由于采用了以上技术方案，本发明的优点在于：

[0019] (1)充分利用了小波变换在时间、频率分辨率上均可调整的局部化时频分析的能力，利用离散二进小波变换系数的模极大值的位置同信号的局部奇异性密切相关性，辨识出信号的局部奇异点，可以定量分析出异常地质体的边界信息，减少了人为经验的影响，更为客观。

[0020] (2)直方图均衡化的方法能提高图像的局部差异对比，增强细节显示效果，与小波变换方法结合，能够使图像信息更为丰富和直观，提高了不良地质体的边界分辨精度。

[0021] 本发明的基于二维视电阻率数据利用小波变换和图像增强表达相结合的数据挖掘方法进行地质异常体识别的方法，能够对物探数据进行有效分析，获得对断层、破碎带、溶洞等异常区域更为精细化的描述，为设计、施工提供科学依据。

[0022] 附图表说明

[0023] 附图1 2号测线上的部分二维视电阻率数据原始灰度图。

[0024] 附图2附图1区段进行直方图均衡化后的结果，其细节对比度明显增强。

[0025] 附图3高阶异常分界点回溯定位。

[0026] 附图4异常分界点与直方图均衡化后的二维视电阻率数据灰度图叠置效果图(小波阶 数为4)。

[0027] 附图5异常分界点与直方图均衡化后的二维视电阻率数据灰度图叠置效果图(小波阶数为3)。

具体实施方式

[0028] 下面结合附图及实施例对本发明的利用二维视电阻率数据识别地质异常体方法进一步详细说明，见附图。

[0029] 一种利用二维视电阻率数据识别地质异常体的方法，基于二维视电阻率数据利用小波变换和图像增强表达相结合的数据挖掘方法进行地质异常体识别，其识别方法按以下步骤进行：

[0030] ①利用灰度变换函数T(r)将研究区域内用物探方法获得的任一测线上所有测点的二维视电阻率数据数组a[x][y]进行直方图均衡化，均衡化的二维视电阻率数据为s[x][y]。

[0031] 其中，x为测点桩号，y为该测点下某一深度范围内任一点的深度。

[0032] ②利用小波变换检测二维视电阻率数据的奇异性：

[0033] a)将任意桩号位置xp随深度变化的视电阻率数据a[xp][y]看成是随着深度y变化的一维信号f(y)；

[0034] 其中，p＝1,2,3…m。

[0035] b)设f(y)为输入信号，θ(y)是某一起平滑作用的低通平滑函数，ψ(y)为θ(y)的一阶导数 记 则对尺度为2j的小波变换为：

[0036]

[0037] 其中，j为小波变换的阶数。

[0038] c)利用Matlab软件的Findpeaks函数找出小波变换 的模极大值，确定第j阶的模极大值点的位置序列

[0039] 其中，设第j阶有n个极大值点，k＝1,2,3…n。

[0040] d)对第j阶的每个极大值点位置 向下在相邻尺度j-1阶上寻找所对应的极大值传播点位置 然后再以j-1阶的 为模极大值点的位置序列再次向下寻找j-2阶上所对应的极大值传播点位置 继续向下回溯，直到寻找到第一阶上的 为止。得到的第一阶上  的位置即为第j阶的模极大值点的位置序列 所对应的异常分界点所在的位置。

[0041] e)对其他所有桩号的视电阻率数据重复b)～d)步骤执行，直到所有桩号视电阻率数据均分析完成。

[0042] ③将①中直方图均衡化后的视电阻率数据图像与②中利用小波变换检测出的奇异点叠 加显示，如果奇异点成团、成线状分布，并且视电阻率图像也呈明显低阻特征，可据此勾勒出异常地质体边界。

[0043] 经过反复实验发现，所述的小波变换的阶数j＝4时，识别出的奇异点大都能与视电阻率图像的低阻部位相对应，而且地质资料中断层带、钻孔涌水等异常部位都能在该方法获得的异常区域上有所反映，识别地质异常体的效果更好。

[0044] 所述的小波变换的阶数j＝3时，相较4阶结果，识别出的奇异点明显增多，二维视电阻率数据奇异点检测的敏感程度更高。

[0045] 异常地质体边界实际上就是视电阻率曲线中的剧烈变化之处，检测异常边界点实质上就是检测视电阻率数据中的奇异点和不规则的突变部分。小波变换是时频分析的一种工具，被称为数学分析上的“放大镜”和“显微镜”。不但具备有局部化时频分析的能力，而且时间分辨率和频率分辨率均可以调整，当前信号突变检测已经成为小波变换应用的一个重要方面。离散二进小波变换系数的模极大值的位置同信号的局部奇异性密切相关，利用这些极大值可以辨识出信号的局部奇异点，本发明正是利用小波变换的这种特性来检测视电阻率数据中的异常分界点。

[0046] 岩石中的破碎带、裂隙、溶蚀、岩性变化、地下水等会引起视电阻率的急剧下降，特别是在与地下水的相互作用下。本发明采用图像增强表达的方法将视电阻率数据的高低在灰度图上直观显示，再与利用小波变换检测的异常分界点数据叠置来综合识别异常地质体。

[0047] 具体实施例

[0048] 下面结合实施例对本发明的利用二维视电阻率数据识别地质异常体方法进一步详细介绍：

[0049] 实施例一

[0050] 石鼓望城坡香炉山隧洞是滇中引水最关键的控制性工程。该隧洞距离长，埋深大，在施工过程中极有可能发生突泥突水等问题，对施工安全造成巨大的隐患。为了分析前期勘探所获得的大地电磁物探数据(二维视电阻率数据)中的不良地质体信息，采用一种利用二维视电阻率数据识别地质异常体方法，它基于二维视电阻率数据利用小波变换和图像增强表达相结合的数据挖掘方法进行地质异常体识别，其识别方法按以下步骤进行：

[0051] ①利用灰度变换函数T(r)将研究区域内用物探方法获得的任一测线上所有测点的二维视电阻率数据数组a[x][y]进行直方图均衡化，其中，x为测点桩号，y为该测点下某一深度范围内任一点的深度。实例数据选择2号测线桩号从DL Ⅱ 000+000至DL Ⅱ 055+960区段的视电阻率数据，全区段长55960m，在测线上每隔10m取一个测点桩号来进行分析，深度上从地表开始，以10m为间距取点，最深深度达1420m，形成的二维视电阻率数据数 组部分见附表1所示，其部分原始灰度图如附图1所示。利用matlab软件的histeq函数对视电阻率数据进行直方图均衡化，附图2为附图1区段进行直方图均衡化后的结果，其细节对比度明显增强。

[0052] 附表1

[0053]

[0054]

[0055]

[0056]

[0057] ②利用一维二进小波变换检测二维视电阻率数据的奇异性：

[0058] a)将任意桩号位置xp随深度变化的视电阻率数据a[xp][y]看成是随着深度y变化的一维信号f(y)；

[0059] 其中，p＝1,2,3,…,5596。

[0060] b)设f(y)为输入信号，选用高斯函数 作为平滑函数，尺度因子a＝2j，ψ(y)为θ(y)的一阶导数 记 则对尺度为2j的小波变换为：

[0061]

[0062] 其中，j为小波变换的阶数。以DL Ⅱ 049+500测点桩号为例，这里最大阶数j取4，计算出的各阶小波变换值见附表2所示。由于求小波变换的过程中要求一阶导数，因此，高程1810处无小波值。

[0063] 附表2

[0064]高程(m) 一阶小波值 二阶小波值 三阶小波值 四阶小波值
2680 666.53 693.04 726.19 715.42
2670 831.43 820.98 817.42 768.12
2660 991.33 948.50 907.92 818.47
2650 1139.00 1071.15 995.54 865.55
2640 1269.79 1185.13 1078.05 908.45

[0065]2630 1380.92 1287.50 1153.31 946.30
2620 1471.03 1376.26 1219.26 978.30
2610 1540.29 1450.27 1274.01 1003.75
2600 1590.81 1508.92 1315.89 1022.05
2590 1626.52 1551.76 1343.47 1032.78
2580 1652.12 1578.07 1355.63 1035.65
2570 1670.98 1586.55 1351.61 1030.54
2560 1682.88 1575.19 1331.04 1017.55
2550 1682.20 1541.50 1294.03 996.91
2540 1658.05 1482.99 1241.20 969.07
2530 1596.64 1397.99 1173.72 934.63
2520 1485.89 1286.44 1093.31 894.34
2510 1320.61 1150.56 1002.20 849.08
2500 1106.24 995.14 903.08 799.84
2490 859.32 827.32 798.96 747.68
2480 603.86 655.82 693.04 693.70
2470 365.08 489.84 588.52 639.01
2460 163.06 337.88 488.43 584.69
2450 9.09 206.68 395.47 531.79
2440 -94.62 100.51 311.84 481.25
2430 -152.76 20.96 239.19 433.94
2420 -174.04 -32.81 178.55 390.57
2410 -168.54 -63.53 130.34 351.75
2400 -145.73 -75.10 94.44 317.94
2390 -113.36 -71.86 70.24 289.45
2380 -77.43 -57.95 56.83 266.45
2370 -42.30 -36.88 53.04 248.98
2360 -10.54 -11.34 57.61 236.98
2350 17.44 16.90 69.25 230.25

[0066]2340 43.25 46.81 86.75 228.53
2330 69.55 77.94 109.00 231.46
2320 98.88 110.18 135.05 238.63
2310 132.79 143.57 164.06 249.60
2300 171.15 178.19 195.33 263.91
2290 211.61 214.14 228.26 281.06
2280 250.36 251.56 262.33 300.57
2270 284.56 290.54 297.08 321.96
2260 315.50 331.07 332.07 344.77
2250 349.48 372.85 366.91 368.56
2240 394.54 415.09 401.25 392.91
2230 454.52 456.45 434.77 417.42
2220 524.30 495.22 467.18 441.74
2210 590.22 529.70 498.28 465.52
2200 635.96 558.77 527.90 488.44
2190 650.89 582.36 555.89 510.22
2180 636.36 601.70 582.14 530.57
2170 606.12 619.08 606.53 549.23
2160 580.78 637.28 628.87 565.93
2150 578.89 658.59 648.96 580.43
2140 609.58 684.03 666.50 592.45
2130 670.07 712.78 681.16 601.70
2120 748.19 742.11 692.62 607.88
2110 827.52 767.88 700.60 610.62
2100 891.82 785.52 704.95 609.54
2090 927.55 791.07 705.68 604.16
2080 925.07 782.28 702.98 593.97
2070 880.34 759.31 697.27 578.35
2060 797.26 724.84 689.07 556.64

[0067]2050 689.43 683.79 678.93 528.09
2040 578.57 642.31 667.27 491.91
2030 488.84 606.74 654.17 447.26
2020 438.62 582.38 639.18 393.31
2010 434.20 572.62 621.08 329.28
2000 468.76 578.49 597.75 254.46
1990 527.26 598.54 566.05 168.28
1980 594.28 629.00 521.85 70.42
1970 660.61 663.89 460.13 -39.22
1960 725.33 694.89 375.34 -160.36
1950 792.47 710.94 261.76 -292.36
1940 863.72 697.71 114.14 -434.16
1930 930.49 637.35 -71.61 -584.21
1920 969.04 509.01 -297.68 -740.54
1910 941.01 290.53 -563.64 -900.74
1900 798.60 -38.21 -865.89 -1062.03
1890 490.22 -490.55 -1197.29 -1221.31
1880 -36.49 -1067.35 -1547.20 -1375.28
1870 -831.17 -1751.65 -1901.90 -1520.55
1860 -1924.99 -2505.35 -2245.39 -1653.78
1850 -3293.22 -3269.88 -2560.57 -1771.83
1840 -4807.13 -3972.02 -2830.73 -1871.86
1830 -6214.37 -4534.84 -3041.03 -1951.51
1820 -7190.67 -4891.77 -3179.98 -2008.95

[0068] c)利用Matlab软件的Findpeaks函数找出小波变换 的模极大值，确定出第j阶的模极大值点的位置序列 其中，设第j阶有n个极大值点，k＝1,2,3…n。实例中桩号为DLⅡ049+500测点处视电阻率数据经小波变换后的各阶模极大值的位置列表如附表3所示。

[0069] 附表3

[0070]小波阶数 模极大值的位置(m)
一阶小波 1920  2010  2090  2150  2190  2420  2560
二阶小波 1950  2010  2090  2400  2570
三阶小波 2090  2370  2580
四阶小波 2110  2340  2580

[0071] d)对第j阶的每个极大值点位置 向下在相邻尺度j-1阶上寻找所对应的极大值传播点位置 然后再以j-1阶的 为模极大值点的位置序列再次向下寻找j-2阶上所对应的极大值传播点位置 继续向下回溯，直到寻找到第一阶上的 为止。得到的第一阶上  的位置即为第j阶的模极大值点的位置序列 所对应的异常分界点所在的位置。桩号为DL Ⅱ 049+500测点处，第4阶共有3个模极大值点，分别在高程2580m、2340m和2110m处，对每个极大值点向下在相邻第3阶上寻找到所对应的极大值传播点，分别为高程2580m、2370m和2090m处，再向二阶回溯，得到对应的异常点位置为高程2570m和2090m处，再次回溯到一阶，得到对应的异常点位置为高程2570m、2400m和2090m处，再次回溯到一阶，得到异常点位置为高程2560m、2420m和2090m处，回溯定位过程如附图3所示。

[0072] e)对其所有他桩号的纵向视电阻率数据按b)～d)步骤执行，直到所有桩号视电阻率数据均分析完成。

[0073] ③将①中直方图均衡化后的视电阻率数据图像与②中利用小波变换检测出的奇异点叠加显示，如果奇异点成团、成线状分布，并且视电阻率图像也呈明显低阻特征，可据此勾勒出异常地质体边界。附图4为测点桩号为DL Ⅱ 049+300至DL Ⅱ 052+500区段，利用上述方法，当小波阶数j＝4时处理得到的结果；附图5为相应区段在j＝3时处理得到的结果。

[0074] 所述的小波变换的阶数j＝4时，识别出的奇异点大都能与视电阻率图像的低阻部位相对应，而且地质资料中断层带、钻孔涌水等异常部位都能在该方法获得的异常区域上有所反映，识别地质异常体的效果更好。而所述的小波变换的阶数j＝3时，相较4阶结果，识别出的奇异点明显增多，二维视电阻率数据奇异点检测的敏感程度增高。

[0075] 从异常界面点的展布看，大部分段落中都存在与地层起伏基本一致、几近连续可以连成线的异常界面点串，基本上可以反映出地层的风化界面，上覆地层由于卸荷、地表水、风化侵蚀作用，导致视电阻率值较低；异常界面线下的地层相对来视电阻率值较高，岩体较为完整。