2017, Vol. 60
可控源音频大地电磁法(Controlled Source Audio-frequency Magnetotelluric Method，CSAMT)是在大地电磁法和音频大地电磁法的基础上，针对大地电磁法场源的随机性强和信号弱的特点而提出来的一种人工源频率域测深方法，具有信噪比高、探测深度大、空间分辨率高、工作效率高等优点.CSAMT方法在地热资源探测以及矿产资源和地下水资源勘查等方面，都取得了良好的效果(李章喜等，2006；李冰等，2013；于昌明，1998；何俊飞，2013；冀显坤等，2014；祝杰，2012).但是在CSAMT信号采集过程中，地表或近地表往往存在局部导电不均匀体，当电流流过不均匀体的表面时，就会形成电荷累积效应，导致电场发生畸变，产生静态效应.在“视电阻率-频率”双对数坐标系中相邻各测点的视电阻率曲线沿着视电阻率轴上下移动，反映在视电阻率拟断面图上，会出现横向范围不大的陡立状密集带，掩盖地下由矿体或目标物引起的异常，影响解释精度，需要进行静态校正.
静态校正一般有以下几种常用方法(黄兆辉等，2006)：曲线平移法、空间滤波法、联合反演法、电磁阵列剖面(EMAP)方法、数值模拟方法、小波滤波方法等.曲线平移法(黄潜生等，2004)比较依赖于处理人员识别曲线的能力；空间滤波法(李金铭和罗延钟，1996；杨生，1994；阎述和陈明生，1996)校正的关键在于滤波窗口的宽度和滤波系数的选择；联合反演法(张旭等，2010)相当于重复做了两种地球物理方法，成本高，实用性不强.小波变换(宋守根等，1995)是一种变窗口变换，将地表不均匀性引起的静态效应与大构造异常在数学上进行直观的区分，有较高的分辨率和校正精度.
宋守根等(1995) 首次利用小波分析方法来进行电磁测深中静态效应的识别、分离与压制，为用小波分析研究静态校正奠定了理论基础；陈辉(2007) 以及程云涛(2008) 先后将小波方法应用到CSAMT静态校正上.这些研究都是将视电阻率数据进行多尺度分解后的所有细节部分的小波系数设置为零，因此会损失一部分大构造异常的信息，使重构后的信号发生畸变.由于在小波变换的过程中，静态效应对小尺度上小波系数的影响比较大，而在大尺度上，小波系数主要由大构造异常控制.因此本文提出一种将小波变换模极大值法与阈值法相结合的静态效应压制方法，该方法将视电阻率数据作多尺度分解后，先求取不同尺度的模极大值，然后设置相应的阈值来压制静态效应.相比于将所有细节部分设置为零的方法，本文提出的方法在压制静态效应的基础上，还可以最大限度保留大构造异常信息.
2 静态效应的识别 2.1 小波变换及模极大值连续小波变换的定义如式(1) 所示：
|
(1) |
式中ψa，τ(t)为小波基函数，a为尺度因子，τ为平移因子，变换结果Wf(a，τ)称为小波变换系数.Wf(a，τ)可以视为信号f(t)在基矢量ψa，τ(t)方向上的投影的大小，当基矢量不断变化时，小波系数就相当于原信号在不同基矢量方向上分量的大小(刘涛等，2006)，小波变换示意图如图 1所示.
|
图 1 小波变换示意图 Fig. 1 Sketch of wavelet transformation |
若对x0邻域内的任意一点，有Wf(s0，x0)≥Wf(s0，x)，则称(s0，x0)为f(x)在尺度s0上的小波变换模极大值点，Wf(s0，x0)为f(x)在尺度s0上的小波变换模极大值.
2.2 李氏指数若信号在某一点或某一个区间内可微的次数越高，那么该信号越平滑，奇异性越弱.在数学中，用Lipschitz指数(以下简称李氏指数)来定量描述函数的奇异性.信号x(t)在t0点的邻域(t0-h，t0+h)内的泰勒级数(胡广书，2004)是:
|
(2) |
Tt0(t)和x(t)的误差et0(t)满足下列关系:
|
(3) |
式(3) 给出了误差函数的上界，其中n为整数，而李氏指数利用一个非整数α来描述这个误差的上界，即若存在常数k>0以及多项式Tt0(t)，使得：
|
(4) |
称x(t)在t0处具有李氏指数α.
函数f(x)在区间(a，b)上具有一致李氏指数α的充要条件为：存在常数k>0，使得对所有的x∈(a，b)，f(x)的小波变换Wf(s，x)满足不等式(5) (Mallat and Hwang，1992)：
|
(5) |
其中k为常数，s为尺度. 当进行二进离散小波变换时，式(5) 变为：
|
(6) |
对式(6) 两边取对数得到：
|
(7) |
对于任意两个相邻尺度si=2j，si+1=2j+1分别代入式(7) 并相减得：
|
(8) |
式(8) 中x1是x0在尺度j时的传播点.
寻找小波变换模极大值点在不同尺度间传播点的步骤如下(张玉新等，2009):
(1) 设x0是尺度J上的模极大值点，x1、x2是x0前后相邻的2个模极大值点，x′1是x1传播到J-1尺度上的对应模极大值点，则x0在J-1尺度上对应的传播点将在区间[max(x1，x′1)，x2]搜索；
(2) 在[max(x1，x′1)，x2]区域内，若存在x′0=x0，且满足W2J-1f(x′0)和W2Jf(x0)符号相同，则x′0是x0的传播点；若不存在这样的点，则在该尺度寻找与W2Jf(x0)符号相同且最接近的那个模极大值点作为x′0，即满足式(9) ：
|
(9) |
(3) 若在区间[max(x1，x′1)，x2]内没有找到符合以上条件的对应传播点，根据Mallat理论，则将x0点予以剔除；
(4) 重复以上步骤，直到所需的尺度.
由于大构造异常的小波变换模极大值随尺度的增加而增加，由公式(8) 可以看出李氏指数一般为正值，而由静态效应引起的小波变换模极大值随尺度的增加而减少，由公式(8) 可知李氏指数为负值，所以可以根据计算所得的李氏指数判断奇异点是否为静态效应，根据第J尺度上的小波变换模极大值点在第一尺度上的传播点判断出现静态效应的位置.
3 静态效应的压制 3.1 小波分解对视电阻率数据进行二进离散小波变换时，小波基和分解尺度的选择会影响校正效果.对于具体的被测信号，要选择合适的小波母函数.一般来讲，对于检测有阶跃变化的信号易采用反对称小波，而对于检测具有尖峰脉冲变化的信号宜采用对称小波(刁彦华等，2004).
分解尺度过大，将会失去有用信号的一些重要特征，分解尺度过小，静态效应的模极大值衰减不充分，在实际操作中，可以令
随着尺度的增加，大构造异常控制的模极大值点占据优势，而静态效应的模极大值点主要分布在小尺度上，所以随着尺度的增加，应适当地减小阈值，这样既可以压制静态效应，又能很好地保留大构造异常信号.
参照彭园园(2011) 的文献，并结合CSAMT静态校正存在的问题，确定本文阈值设置策略如下：在最大尺度上阈值设置为thr=c×M/J，在分解尺度1~J－1上设置阈值thr=c×M2/ln(j+1) ，其中c一般取0.8，M为相对应尺度上所有模极大值中的最大值.
阈值设置完成后，对于不同尺度上的小波变化模极大值采取硬阈值方式处理，即对大于等于该尺度阈值的小波变换模极大值予以保留，小于阈值的小波变换模极大值设置为零.
3.3 Mallat交替投影法重构如果原始信号f(x)存在于两个凸集D1、D2的交集Λ中，即Λ=D1∩D2，f(x)∈Λ.那么对小波系数的重构即可化为至少求出Λ中的一个元素.
如果投影到两个凸集D1、D2上的投影算子分别是P1、P2，那么定义一个变量P：P=P1·P2，重构信号就是通过多次重复投影P来近似逼近原始信号f(x)，f(x)的估计值
|
(10) |
式(10) 中f(x0)为凸集D1中的最初估计.重复过程如图 2所示.因为二进离散小波变换所需的数据量是按照2的指数倍增长的，所以当测点数比较少时，在小波分解前需要做一步插值处理，基于小波变换模极大值法与阈值法的静态校正流程图如图 3所示.
|
图 2 交替投影法示意图 Fig. 2 Schematic diagram of the alternating projection method |
|
图 3 静态校正流程图 Fig. 3 Flowchart of static correction |
正演模型一如图 4所示，在测点20~22之间有一个顶部埋深为零，尺寸大小为100 m×100 m的静态体，电阻率为500 Ωm.距地表 900 m的范围内为均匀大地，电阻率为200 Ωm.在距地表 600 m处有一个电阻率为1000 Ωm的高阻体，宽度为400 m，高度为300 m，测点间距为50 m.
|
图 4 正演模型一示意图 Fig. 4 Schematic diagram of forward model I |
图 5为正演数值模拟ρa拟断面图，从图 5可以明显看出在测点20~22之间有陡立状的等值线出现，经计算在测点20和22左右李氏指数近似值为-0.28，表明了静态体的存在.图 6为将所有尺度细节部分的小波系数设置为零后的ρa拟断面图，从图 6可以看出，静态效应得到了一定的压制，但是由于过度平滑，导致大构造异常体的横向范围严重变宽，与实际的地质情况差别很大.图 7为利用本文提出的模极大值和阈值结合方法处理后的ρa拟断面图，与图 6相比，不仅静态效应压制效果更好一点，并且很好地保留了大构造异常体的信息，与模型基本相符.
|
图 5 模型一正演数值模拟ρa拟断面图 Fig. 5 ρa pseudo-section from forward numerical simulation on model I |
|
图 6 模型一所有细节系数设置为零后的小波滤波结果 Fig. 6 Wavelet filtering result after all detail coefficients are set to zero on model I |
|
图 7 模型一对小波变换模极大值进行阈值处理后的结果 Fig. 7 The result after the wavelet transformation modulus maxima are processed by thresholds on model I |
正演模型二如图 8所示，该模型设置20个测点，点距为50 m，在地下0~50 m范围内有一盖层，电阻率为100 Ωm，在测点10附近有一个顶部埋深为50 m，尺寸大小为50 m×50 m的低阻静态体，电阻率为10 Ωm，在测点16附近有一个顶部埋深为50 m，尺寸大小为50 m×50 m的高阻静态体，电阻率为500 Ωm，距离地面50~950 m范围内的背景电阻率为50 Ωm，然后在距地面650 m处测点14~18之间有个宽度为200 m、厚度为300 m的高阻体，电阻率为100 Ωm.
|
图 8 正演模型二示意图 Fig. 8 Schematic diagram of forward model II |
图 9为正演数值模拟ρa拟断面图，从图 9可以看出，由于两个静态体的存在，导致测点10附近出现陡立状的等值线，而测点16附近的盖层和下方高阻体连为一体，无法辨别.经计算，该两点的李氏指数分别为-0.32和-1.86.图 10为将所有尺度细节部分的小波系数设置为零后的ρa拟断面图，从图 10可以看出，静态效应得到了一定的压制，但是测点1~13范围内的电阻率值过度降低，没有呈现出均匀的空间现象.图 11为利用本文提出的模极大值和阈值结合方法处理后的ρa拟断面图，与图 10相比，不仅静态效应压制效果更好一点，并且很好地保留了大构造异常体的信息，将模型二的信息很好地体现出来.
|
图 9 模型二正演数值模拟ρa拟断面图 Fig. 9 ρa pseudo-section from forward numerical simulation on model II |
|
图 10 模型二所有细节系数设置为零后的小波滤波结果 Fig. 10 The wavelet filtering result after all detail coefficients are set to zero on model II |
|
图 11 模型二对小波变换模极大值进行阈值处理后的结果 Fig. 11 The result after the wavelet transformation modulus maximum are processed by thresholds on model II |
2013年5月22日起“地面电磁探测(SEP)系统研制——野外试验研究”项目组在内蒙古兴和县曹四夭钼矿区开展野外矿区的实际试验应用.
5.1 矿区地球物理特征及工作区采集情况曹四夭矿区的激电异常呈南北向条带状和串珠状展布，并且存在有多个高异常值中心，其视幅频率(Fs值)均大于1.5%，最高值为14.75%，视电阻率值为50~400 Ωm，最高值为500 Ωm(李香资等，2012).矿区范围内，激电异常的高极化率和低电阻率值对应得很好，暗示了深部含硫化物地质体的存在.曹四夭矿区的激电异常、钨-钼-砷-铜-锌异常与含钨钼热液蚀变带相互叠加的部位是找矿的有利地段.
根据已有地质资料及课题组要求，采用V8系统进行了勘探，工作区布置6条勘探线，分别呈南北向、东西向及北西向各2条，本文选取东西向H08线勘探数据为例，进行静态校正.H08线剖面长度3.5 km，点距25 m，线上有140个电测深点.工作方式每个排列5个电道和1个磁道，发射极距1 km，收发距10 km，发射电流16 A.
5.2 实测资料处理结果图 12为由H08测线原始数据绘制的视电阻率断面图.由图 12可以看出，在测点90、120、130附近存在严重的静态效应.经计算，这三个测点附近的李氏指数近似值分别为-2.32、-0.22、-1，表明了静态效应的存在.
|
图 12 H08测线视电阻率断面图 Fig. 12 Apparent resistivity section of line H08 |
图 13为H08原始数据的反演结果，由图 13可以看出，由于静态效应的存在，导致在测线的2000~2500 m、3000~3500 m范围内出现陡立的等值线，呈现虚假的高阻，掩盖了实际的地质情况.针对H08测线的静态效应，采取小波变换模极大值和阈值相结合的方法压制静态效应.首先，将140个测点数据采用三次样条插值方法增加到512个测点数据，然后利用二进离散小波变换对数据进行八尺度分解，求取每一尺度的小波变换模极大值，计算相应尺度阈值，然后采用硬阈值方法处理小波变换模极大值，最后利用Mallat交替投影法进行重构.图 14为静态校正后的反演结果.由图 14可以看出，静态效应得到消除，测线1200~2500 m范围内的矿体(电阻率值为50~500 Ωm，图中绿色部分)明显，与实际钻井资料一致.
|
图 13 H08测线原始数据反演结果 Fig. 13 The inversion result using original data of line H08 |
|
图 14 H08测线静态校正后反演结果 Fig. 14 The inversion result using data after static correction of line H08 |
小波分析方法在CSAMT静态效应校正中起到了很好的作用，但在压制的过程中，也会引起校正不足或者校正过度的问题.本文在小波分析的基础上，提出了基于小波变换模极大值和阈值法相结合的方法进行静态效应校正.首先，对获得的测点-视电阻率数据进行多尺度小波变换，求取每一尺度上的小波变换模极大值，进一步计算李氏指数来判断静态效应.对于静态效应的压制，区别于将所有细节系数都设置为零的做法，本文在每一尺度上通过设置合理的阈值，对相应尺度的小波变换模极大值利用硬阈值函数进行处理，然后用Mallat交替投影法进行重构即可获得静态压制后的视电阻率数据.与将所有的细节系数设置为零相比，小波变换模极大值和阈值结合策略在压制静态效应的同时，可以很好地保留大构造异常的信息.通过正演模型和实际案例验证了本文方法的有效性.
致谢感谢吉林大学地球信息探测仪器教育部重点实验室提供了研究平台，感谢吉林大学仪器科学与电气工程学院电磁组的全体成员给予的帮助和指导.
| Chen H. 2007. The research on the technology and method of the static effect of controlled source audio frequency magnetotelluric method[Master's thesis] (in Chinese). Beijing:China University of Geosciences. | |
| Cheng Y T. 2008. The identification and correction of CSAMT static effect[Master's thesis] (in Chinese). Changsha:Central South University. | |
| Diao Y H, Wang Y T, Chen G T. 2004. Singularity detection of signals based on wavelet transform modulus maximum. Hebei Journal of Industrial Science & Technology (in Chinese), 21(1): 1-3. | |
| He J F. 2013. The application of controlled source audio-frequency magnetotelluric method in Fankou Lead-zinc deposit. Chinese Journal of Engineering Geophysics (in Chinese), 10(6): 763-770. | |
| Hu G S. Modern Signal Processing Tutorial (in Chinese).Beijing: Tsinghua University Press, 2004: 372-373. | |
| Huang Q S, Wang Y S, Wang W M. 2004. The application of magnetotelluric curve correction technique in Six Panshan Basin. Journal of Jianghan Petroleum Institute (in Chinese), 26(3): 76-78. | |
| Huang Z H, Di Q Y, Hou S L. 2006. CSAMT static correction and its application. Progress in Geophysics (in Chinese), 21(4): 1290-1295. | |
| Ji X K, Wang X D, Tang S S, et al. 2014. The application of CSAMT and symmetric quadrupole sounding method to groundwater survey. Chinese Journal of Engineering Geophysics (in Chinese), 11(3): 342-345. | |
| Li B, Wang X Q, Ding Y H, et al. 2013. Application study of CSAMT method to investigation and evaluation project of bauxite. Journal of Engineering Geology (in Chinese), 21(6): 987-994. | |
| Li J M, Luo Y Z. New Advances in Electrical Prospecting (in Chinese).Beijing: Geological Publishing House, 1996: 37-41. | |
| Li X Z, Ban Y H, Quan Z X, et al. 2012. Discuss on the molybdenum deposit geochemical characteristics and metallogenic model in Xinghe County, Inner Mongolia. Geological Survey and Research (in Chinese), 35(1): 39-46. | |
| Li Z X, Gu H M, Ran C G. 2006. Application of CSAMT method to geological exploration in the Motian mountain tunnel. Chinese Journal of Engineering Geophysics (in Chinese), 3(3): 187-191. | |
| Liu T, Zeng X L, Zeng J. Introduction to Practical Wavelet Analysis (in Chinese).Beijing: National Defense Industry Press, 2006: 40-45. | |
| Mallat S, Hwang W L. 1992. Singularity detection and processing with wavelets. IEEE Transactions on Information Theory, 38(2): 617-643. DOI:10.1109/18.119727 | |
| Peng Y Y. 2011. The application of wavelet analysis in one-dimensional signal de-noising[Master's thesis] (in Chinese). Beijing:University of Posts and Telecommunications. | |
| Song S G, Tang J T, He J S. 1995. Wavelets analysis and the recognition, separation and removal of the static shift in electromagnetic soundings. Chinese Journal of Geophysics (Acta Geophysica Sinica) (in Chinese), 38(1): 120-128. | |
| Yan S, Chen M S. 1996. Static offset and correction in frequency domain electromagnetic sounding. Oil Geophysical Prospecting (in Chinese), 31(2): 238-247. | |
| Yang S. 1994. The EMAP data processing technique and its application to the CSAMT survey. Geological Exploration for Non-ferrous Metals (in Chinese), 3(6): 345-348. | |
| Yu C M. 1998. The application of CSAMT method in looking for hidden gold mine. Chinese Journal of Geophysics (Acta Geophysica Sinica) (in Chinese), 41(1): 133-138. | |
| Zhang X, Liu K H, Li X. 2010. Application of static correction CSAMT-joint inversion method. Northwestern Geology (in Chinese), 43(2): 38-43. | |
| Zhang Y X, Teng G F, Zhao Y, et al. 2009. Research on wavelet modulus maxima in the signal de-noising method. Journal of Agricultural University of Hebei (in Chinese), 32(1): 114-116. | |
| Zhu J. 2012. The applied technology of controlled source audio magnetotelluric on geothermal field[Master's thesis] (in Chinese). Chengdu:Chengdu University of Technology. | |
| 陈辉. 2007. 可控源音频大地电磁法静态效应校正技术及其方法研究[硕士论文]. 北京:中国地质大学. | |
| 程云涛. 2008. CSAMT静态效应的识别及校正[硕士论文]. 长沙:中南大学. | |
| 刁彦华, 王玉田, 陈国通. 2004. 基于小波变换模极大值的信号奇异性检测. 河北工业科技, 21(1): 1–3. | |
| 何俊飞. 2013. 可控源音频大地电磁法在凡口铅锌矿中的应用. 工程地球物理报, 10(6): 763–770. | |
| 胡广书. 现代信号处理教程.北京: 清华大学出版社, 2004: 372-373. | |
| 黄潜生, 王友胜, 汪卫毛. 2004. 大地电磁曲线校正技术在六盘山盆地研究中的应用. 江汉石油学院学报, 26(3): 76–78. | |
| 黄兆辉, 底青云, 侯胜利. 2006. CSAMT的静态效应校正及应用. 地球物理学进展, 21(4): 1290–1295. | |
| 冀显坤, 王小多, 唐圣松, 等. 2014. CSAMT法和对称四极测深法在地下水调查中的应用. 工程地球物理学报, 11(3): 342–345. | |
| 李冰, 王秀全, 丁云河, 等. 2013. CSAMT法在铝土矿调查评价项目中的应用研究. 工程地质学报, 21(6): 987–994. | |
| 李金铭, 罗延钟. 电法勘探新进展.北京: 地质出版社, 1996: 37-41. | |
| 李香资, 班宜红, 权知心, 等. 2012. 内蒙古兴和县曹四夭钼矿床地球化学特征及成矿模型探讨. 地质调查与研究, 35(1): 39–46. | |
| 李章喜, 顾汉明, 冉昌国. 2006. CSAMT在摩天岭隧道地质勘察中的应用. 工程地球物理学报, 3(3): 187–191. | |
| 刘涛, 曾祥利, 曾军. 2006. 实用小波分析入门. 北京:国防工业出版社, 40. | |
| 彭园园. 2011. 小波分析在一维信号去噪中的应用[硕士论文]. 北京:北京邮电大学. | |
| 宋守根, 汤井田, 何继善. 1995. 小波分析与电磁测深中静态效应的识别、分离及压制. 地球物理学报, 38(1): 120–128. | |
| 阎述, 陈明生. 1996. 频率域电磁测深的静态偏移及校正方法. 石油地球物理勘探, 31(2): 238–247. | |
| 杨生. 1994. EMAP处理方法在CSAMT勘查中的应用. 有色金属矿产与勘查, 3(6): 345–348. | |
| 于昌明. 1998. CSAMT方法在寻找隐伏金矿中的应用. 地球物理学报, 41(1): 133–138. | |
| 张旭, 刘宽厚, 李貅. 2010. CSAMT的静校正应用——联合反演法. 西北地质, 43(2): 38–43. | |
| 张玉新, 滕桂法, 赵洋, 等. 2009. 基于小波变换模极大值的信号去噪方法研究. 河北农业大学学报, 32(1): 114–116. | |
| 祝杰. 2012. 可控源音频大地电磁法在地热勘查中的应用研究[硕士论文]. 成都:成都理工大学. | |