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

Abstract

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

Description

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

技术领域

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

背景技术

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

发明内容

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

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

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

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

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

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

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

b)设f(y)为输入信号，θ(y)是某一起平滑作用的低通平滑函数，ψ(y)为θ(y)的一阶导数 $\mathrm{\ψ}\left(y\right)=\frac{d\mathrm{\θ}}{dy},$ 记 ${\mathrm{\θ}}_{2^j}\left(y\right)=\frac{1}{2^j}\mathrm{\θ}\left(\frac{y}{2^j}\right),$ 则对尺度为2^j的小波变换为：

$W_{2^j}f\left(y\right)=f\left(y\right)*{\mathrm{\ψ}}_{2^j}\left(y\right)=f\left(y\right)*\left(2^j\frac{{\mathrm{d\θ}}_{2^j}\left(y\right)}{dy}\right)=2^j\frac{d}{dy}\[f\left(y\right)*{\mathrm{\θ}}_{2^j}\left(y\right)\]$

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

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

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

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

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

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

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

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

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

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

附图表说明

附图12号测线上的部分二维视电阻率数据原始灰度图。

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

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

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

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

具体实施方式

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

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

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

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

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

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

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

b)设f(y)为输入信号，θ(y)是某一起平滑作用的低通平滑函数，ψ(y)为θ(y)的一阶导数 $\mathrm{\ψ}\left(y\right)=\frac{d\mathrm{\θ}}{dy},$ 记 ${\mathrm{\θ}}_{2^j}\left(y\right)=\frac{1}{2^j}\mathrm{\θ}\left(\frac{y}{2^j}\right),$ 则对尺度为2^j的小波变换为：

$W_{2^j}f\left(y\right)=f\left(y\right)*{\mathrm{\ψ}}_{2^j}\left(y\right)=f\left(y\right)*\left(2^j\frac{{\mathrm{d\θ}}_{2^j}\left(y\right)}{dy}\right)=2^j\frac{d}{dy}\[f\left(y\right)*{\mathrm{\θ}}_{2^j}\left(y\right)\]$

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

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

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

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

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

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

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

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

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

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

具体实施例

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

实施例一

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

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

附表1

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

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

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

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

$W_{2^j}f\left(y\right)=f\left(y\right)*{\mathrm{\ψ}}_{2^j}\left(y\right)=f\left(y\right)*\left(2^j\frac{{\mathrm{d\θ}}_{2^j}\left(y\right)}{dy}\right)=2^j\frac{d}{dy}\[f\left(y\right)*{\mathrm{\θ}}_{2^j}\left(y\right)\]$

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

附表2

高程(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

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.058 -->
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

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.749 -->
					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

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.8610 -->
1830	-6214.37	-4534.84	-3041.03	-1951.51
					1820	-7190.67	-4891.77	-3179.98	-2008.95

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

附表3

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

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所示。

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

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

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

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

Claims (1)

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

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

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

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

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

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

b)设f(y)为输入信号，θ(y)是某一起平滑作用的低通平滑函数，ψ(y)为θ(y)的一阶

导数 $\mathrm{\ψ}\left(y\right)=\frac{d\mathrm{\θ}}{dy},$ 记 ${\mathrm{\θ}}_{2^j}\left(y\right)=\frac{1}{2^j}\mathrm{\θ}\left(\frac{y}{2^j}\right),$ 则对尺度为2j的小波变换为：

$W_{2^j}f\left(y\right)=f\left(y\right)*{\mathrm{\ψ}}_{2^j}\left(y\right)=f\left(y\right)*\left(2^j\frac{{\mathrm{d\θ}}_{2^j}\left(y\right)}{dy}\right)=2^j\frac{d}{dy}\[f\left(y\right)*{\mathrm{\θ}}_{2^j}\left(y\right)\]$

其中，j为小波变换的阶数；

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

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

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

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

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