基于小波分析的重磁数据求导方法及应用.docx

第38卷第6期吉林大学学报(地球科学版) V0138 No62008年1 1月 Journal of Jilin University(Earth Science Edition)Nov2008基于小波分析的重磁数据求导方法及应用方东红1，曾昭发1，陈家林21吉林大学地球探测科学与技术学院，长春 1300262吉林省区域地质矿产调查所。长春 130022摘要：在双正交小波分析下，实现了重磁异常数据高精度的导数计算。该算法是在小波分析得到的系 数基础上，采用MonteCarlo方法进行噪声的估计，根据获得参数分别进行去噪处理。理论计算结果表 明：算法的计算效率高。计算的结果精度高，抗干扰能力强。实际数据处理结果表明：采用基于小渡变换的 方法完成导数处理，结果对地质填图中地质体边界的划分和磁性地层参数分析提供了重要的资料，也为采 用基于导数的位场数据反演解释方法提供了优良的基础数据。关键词：重磁异常数据；导数处理；小波分析；MonteCarlo方法中图分类号：P6312文献标识码：A 文章编号：16715888(2008)061049一06The Derivatives Calculation Based onWavelet Analysis of GM Data and Its ApplicationFANG Dong-hon91，ZENG Zhao-fal，CHEN Jialin2tCollege of GeoExploration Science and TechnologyJilin UniversityChangchun 130026，China2Survey of Regional Geology and Mineral Resources of Jilin Province，Changchun 130022，ChinaAbstract：The wavelet analysis is widely applied in geophysical data processing and interpretation for its multiple resolution analysis，and got very good resultsThe biorthogonal wavelet analysis was used as original basis of wavelet for high accuracy derivates calculation，which serves as important informa tion of anomaly analyses and interpretation，and also as the basic data for further calculation as Euler de convolution，analytic inversion and source parameters interpretation Based on the wavelet coefficients， the Monte-Carlo method was applied to calculate the parameters of noiseThe cases study from real data shows that the algorithm has high accuracy，and its results can help the geological mapping and stratastudyKey words：gravity and magnetic data；derivative processing；wavelet analysis；MonteCarlo method常重要的常规数据处理方法。其结果一方面应用于0引言数据解释和异常分析、研究局部异常和分析地质构重磁数据的求导处理是重磁数据处理和解释非造展布方向 ；另一方面，许多重磁异常解释方法都收稿日期：2008一0318基金项目：国家自然科学基金项目(40774055)作者简介：方东红(1963-)，男，安徽贵池人，博士研究生，主要从事地球物理资料的处理与解释方面的应用研究，Tel：043185159281E-mail：dhfangjlueduca通讯联系人：曾昭发(1966-)，男，广西全州人，教授，主要从事应用地球物理教学与研究，Tel：043188502426。Bmail：zengzfjlueducn。万方数据吉林大学学报(地球科学版)第38卷基于重磁异常的导数值，例如解析信号的解释口 、，(z)=GV2|W，(口，6)似6(z)da万db。(3)磁源参数成像L5，欧拉反褶积【6卅等。莺磁异常导数“的计算精度和质昔商接关系到资料解释和反演结果 式(2)和(3)构成了小波变换的信号分解和重建的变 的精度9。10。由于重磁异常导数处理是一个高通滤 换对。 波处理。其阶数越高，对高频信号的放大倍数越大。12 小波分析求取导数的方法j根据重磁异常的频谱分析，大多数干扰噪声分布在 根据测量的异常数据丁，采用离散小波分析的 高频区域，如果不加控制地对重磁异常进行导数处方法，以双正交小波(daubeehies)构建小波函数，利 理，将直接放大噪声，而有用信号相对受到压制。常用金字塔算法进行分解。根据已有的数学算法 规求导方法是对磁异常进行简单的低通滤波或向上 将位场异常分解为 延拓处理，滤掉高频噪声成分，然后再进行导数处a；。h6；。i。6；。j，+1 6；。i。+，。 理。这种低通滤波在实际应用中很难控制，对高频即AT=C。(4) 有用信号也一起滤掉，影响对高频信号或局部异常 其中：A是小波算子；C是小波算子作用于位场后 的解释。 的小波系数；jmin是根据应用选择小波分解的系数小波变换(wavelet transform)是20世纪90年 所处层；口和b分别是近似系数和细节系数。应用 代发展起来的应用数学分析工具1，它是在傅立叶 对角化算子将C转化为CK：在原始小波基上，进行 变换和Gabor变换基础上发展起来的多尺度的时 导数处理，即KT(在原始异常T上作用一个求导算 频分析工具。由于它具有良好的局部分析功能，广 子K)的系数为：泛应用于图像处理、信号处理等领域1 21。小波分析 CKK。C 7。 (5)的这种多尺度分辨能力也被应用于莺磁数据处理， 式中，K。是K的小波域的算子，其中包括导数处理 并取得了大量的研究结果，主要表现在磁异常滤波、 和去噪处理。其表达式为 微弱信号分析、异常分离等n3。1引。但这些应用往往 云：!|li。=”Ir。(优*c-jmin)，(6) 将分解的系数直接进行取舍，然后进行重构。这种谚=一一15。(7) 取舍只是根据不同的频率成分进行，没有基于异常 式中：仇是一个以27c为周期的函数；K是常数，K一 或噪声的模型，也没有根据模型参数进行处理。本 28。其频域表达式可以选择为 文基于位场的噪声模型，采用MonteCarlo方法计m()一k(一i(扩一1)。(8) 算噪声参数，然后根据这些参数，对不同层次的小波 不同的导数算子K不同：水平导数为石()=(如)。， 系数进行域值处理。 垂向导数为志(甜)=II。；其中，口为求导的阶数。通过小波算法对处理获得的系数进行重构，可1小波分析方法求位场导数方法以得到经过处理的信号：KT=SCK。以上是一维求11 小波变换基本理论1】 导处理的算法。根据这种算子的线性可分性，可以 以一维小波变换为例简单介绍小波分析的基本推广到二维位场异常数据的处理，即石(叫，。)= 理论。定义：若函数妒(z)L1 n L2满足 k(co-)足(眈)，其中叫t，。分别代表两个不同方向的波数17。G一，峄半drooo。(1) I3小波系数的噪声消除方法其中，妒()为妒(z)的傅立叶变换。j小波分析的去噪方式是对小波系数进行域值处 理。即对不同阶次的小波系数，采用不同域值进行令灿(z)=l口Il2妒(孚)，处理。基于两个假设：(1)噪声是正态分布和白噪则函数，(z)L2的小波变换定义为声；(2)分解越细的成分系数中。包含的噪声成分越 W，(口，6)=高。采用MonteCarlo方法估计出不同小波域系l口PJ(z)驴(孚)如。(2) 数的统计参数，并根据这些统计参数确定一系列比 例作为域值的参数。其中，a为尺度因子；6为位置因子。 中心极限定理指出n朝：无论随机变量的分布如 相应的小波逆变换公式为何，它的若干个独立随机变量抽样值之和总是满足万方数据第6期 方东红曾昭发，陈家林：基于小波分析的重磁数据求导方法及应用1051正则分布(即高斯分布)。一个随机变量叩，它满足 分布密度函数厂(z)。将满足分布密度函数厂(z)的 独立随机数相加：Rn=研+】72+巾+ 珈，则R。满 足高斯分布。高斯分布可以由给定的期望值p和方差dr完全确定下来，通常用N(卢t盯2)表示：N(卢，d2)=七exp-(z一)22a2J。(9)oJ 27中心极限定理可以给出蒙特卡洛估计值的偏 差。如果(9)式右边积分的期望值为f，(9)式左边” 次抽样的蒙特卡洛估计值为J。，标准误差为仃，则当 卯充分大时，对任意的A(AO)，有limProb_A磐，。一IA掣)=，呻n住去l1r。e-，佗dt=1_口。 (10)27rJ一1即积分的期望值与蒙特卡洛估计值之差为J。一J lA丛掣。(11)行(11)式的值在误差范围内的概率为口一1，a称图l 基于小波分析方法的位场异常求导方法框图为显著水平。1一a称为置信水平，d为蒙特卡洛估计Fig1 The block diagram of derivative processing of po-值的标准误差，cr2一Vfn。口与A的关系可以由tential anomaly based On wavelet analysis(10)积分公式求得。例如取置信水平1一口=99，可得A=3。这可以解释为：不等式l j。一J I 25)后磁异常曲线见图2a。从图可见，异常的噪3丛圭三成立的概率为99。 声对目标异常影响不大。通过对含噪声异常进行常竹规求导处理，目标体异常的导数响应几乎淹没在噪 由于小波分析的系数含噪声水平不同，采用 声中(图2b)，但采用小波分析方法的目标体导数响MonteCarlo方法对不同频率的系数求取不同的 应非常清楚。统计参数；根据这些不同的统计参数，对高频成分设 为了进一步验证小波分析方法的有效性，将原置相对低的置信水平。而对低频成分设置相对高的 始异常和含噪声异常分别采用小波分析方法进行求置信水平。从而达到消除噪声保留信号的目的。 导处理，结果如图3。从图3可见。含噪声异常通过 图l是该算法的基本框图。 小波分析处理后，两条曲线吻合很好，反映了该方法2 理论计算结果分析对噪声的消除能力和对有效信号的良好保持能力。 采用理论模型进行计算来说明小波分析算法的3 应用分析计算效果。磁场模型由5个垂直磁化的水平圆柱体31 磁异常变化较大的火山岩区域组成，其表达式为(12)53Mc ，=A砉等最端。m， 剖面位于大兴安岭北段的塔河地区1引，地处额尔古纳地块、兴安地块和上黑龙江凹陷交汇部位，地质构造复杂。自新元古代以来，构造一岩浆活动频 A为与磁化率有关的系数，A一2KFS，其中K繁，形成了一系列的侵入岩和火山岩。图4是该剖为水平圆柱体与周围介质磁性差异，F为区域磁场，面的航磁异常图。剖面58 km处强磁异常由中S为圆柱体的截面积。取垂直磁化，并且A=100，酸性黑云母花岗岩引起；基性火山岩英安岩分布在圆柱体中心埋藏深度：h，=36，9，15和20，水平位剖面的1521 km和25-35 km处；其他地方由酸置z。=100，150，200，300和400，剖面采集点数为性火山岩覆盖。由于英安岩的剩磁较强，尽管经过512。原始异常图和叠加商斯噪声(噪声的幅度为万方数据1052吉林大学学报(地球科学版)第38卷bC艇嗽镬点号 点号图2 水平圆柱体理论异常和含噪声异常曲线及其水平一阶导数曲线图Fig2 The curves of magnetic and its noise added of multiple cylinder model and their and horizontal derivative resultsa水平圆柱体模型异常图和含噪声异常曲线lb曲线a的水平一阶导数曲线l墨I小波分析求导方法常规求导办法J似一藏加卜C鑫耗，*睬*翟耗痿啦 趣50 100 150 200 250 300 350 400 450 500、点号图3 原始异常和含噪声异常采用小波方法求导结果比图4 塔河地区航磁异常剖面图较图 Fig4 The profile of magnetic anomalies in Tahe regionFig3 The comparison of derivative results from calcula-ted anomalies and noised magnetic anomalies 噪声，并出现与实际地质情况不相符的假异常。32磁性较弱的沉积岩区域 化极处理，很难能够消除斜磁化对异常解释的影响。为了进一步验证方法的结果，选择大庆地区某图5是采用小波分析方法和常规方法计算的水 剖面进行处理(图6)。剖面地势相对平坦，方向为 平一阶导数、垂直一阶导数和一阶解析异常振幅值。北西，西起巴彦查干乡，东至头台镇南15 km，经过从图可见，采用小波分析方法的计算结果不仅异常了3个重要的构造区。1020 km为中等强磁异常 清晰，而且分辨率高，反映的边界符合野外地质调查 区，经过短暂的过渡，到40100 km处为区域最强的结果。相比较而言，常规方法结果极大地放大了 磁异常区，其后异常逐渐减弱，到负磁异常区。主要加心罄禚心舳O逛蛞O琏霉籁O蹄桩豫O授喜磷静逛塔譬露赫册挺簧嚣 鲫0距离km图5 塔河地区航磁异常水平一阶导数、垂直一阶导数和一阶解析异常振幅 Fig5 The horizontal。vertical derivative and its analytic signal in Tahe region a小漩方法计算结果；b常规方法计算结果万方数据第6期 方东红，曾昭发，陈家林：基于小波分析的重磁数据求导方法及应用1053精度高。 (2)基于小波方法具有较高的计算精度，既消除之了噪声，又有效保留高频信号，为重磁异常的解释和耙酞基于导数的反演解释方法技术提供了基础。弹(3)在两个地区的应用表明，导数处理的结果对地质填图中地质体边界的确定和磁性地层的分析具 有重要的价值。距离km图6 大庆地区航磁异常剖面圈 本文假设噪声为正态分布的白噪声，如果噪声Fig6 The aearomagnetic anomaly profile in Daqing area的模型发生变化，其有效性将受到限制。对不同噪 声模型进行研究，其结果可以进一步服务于本算法。磁性地层是营城组的火山岩地层及下伏的侏罗纪地层。导数处理见图7。从图可见，小波分析结果明 1申宁华，管志宁磁法勘探问题M北京：地质出版社，1985显地反映了区域的磁性变化，与常规方法比较，异常SHEN Ning-hua，GUAN ZhiningThe problems in的稳定性明显提高。magnetic prospectingMBering：Geological Publishing House，1 98521 Nabighian M NToward the three dimensional auto matic interpretation of potential field data via genera- lized Hilbert transforms：fundamental relationsJGeophysics，1984，49：957966IsNabighian M NThe analytic signal of two-dimension al magnetic bodies with polygonal cross-section：its properties and use for automated anomaly interpreta-tionJGeophysics，1972，37：507517距离km204 Nabighian M NAdditional comments on the analytic；signal of two-dimensional magnetic bodies with polygo-15 1毒害釜舅剖器10nal974，39：8592crosssectionJGeophysics，1气蛾罐蛊1553 Smith S，Thurston J B，Da T F i，et a1iSPI“一theimproved sourcemethodJGeo一*0舢parameter imaging*、physical Prospecting，1998，46：14115110。6Nabighian M N，Hansen R OUnification of Euler15and Werner deconvolution in three dimensions via the20406080100120140160180距离kmgeneralized Hilbert transformJGeophysics，2001，66：18051810图7 大庆地区航磁异常采用小波方法和常规方法求取73 Reid A B，Allsop J M，Granser H，et a1Magnetic in垂直、水平一阶导数结果terpretation in three dimensions using Euler deconvolu- Fig7 The vertical and horizontal derivatives of aearomag-tionJGeophysics，1990，55：8091netic anomaly in Daqing area by wavelet method8 Stavrev P YEuler deconvolution using differential and conventional methoda垂直一阶导数lb水平一阶导数similarity transformation占of gravity or magneticanomaliesJGeophysical Prospecting，1 997-45：2072464结论 9谢靖磁源重力异常及其导数的计算方法J长春地质学院学报，1978，8(3)16278(1)基于小波分析的求导方法，配合MenteXIE JingCalculation method of pseudo-gravity anom一Carlo方法进行噪声的估计，算法的效率高，计算的aly and its derivativesJJournal of the Changchun万方数据吉林大学学报(地球科学版) 第38卷Geological Institute，1978，8(3)：6278 油物探，1997，36(2)：8693103曾昭发，吴燕冈，郝立波，等基于泊松定理的重磁异 LI Zong-jieThe application of wavelet transform in 常分析方法及应用J吉林大学学报(地球科学 potential filedDGeophysical Prospecting for Petro 版)。2006，36(2)：279283leum。1997。36(2)：8693ZENG Zhao-fa，wu Yan-gang，HAO Libo，et a115 耿喜哲，刘天佑小波分析在高精度磁测数据处理上The PoissonS theorem based analysis method and ap 的应用I-J石油物探，200039(4)：121124 plication of magnetic and gravity anomaliesJJour-GENG Xizhe，LIU Tian-youThe application of nal of Jilin University(Earth Science Edition)，2006，wavelet analysis in high precious magnetic data pro-36(2)：279283cessingJGeophysical Prospecting for Petroleum。1 13 Mallat S GA theory for muhiresolution signal de2000，39(4)：121124composition：the wavelet representation I-JIEEE 16 Torrence C，Compo GA practical guide to wavelet Trans On Pattern Analysis and Intelligence，1989，IIanalysisJBull Amer Met Soc。1998，79：6178 (7)：67369317 Thomas A。RidsdillSmithThe application of the，12刘贵忠，邸双亮小波分析及其应用M-1西安：西安wavelet transform to the processing of aeromagnetic科技大学出版社。1992，data-DPerths：University of Western Australia， I。lU Gui-Zhong。DI ShuanLiangThe wavelet analy。 2000sis and its applicationMXian：University of D83 Fishman G SMonte Carlo：concepts，algorithms， Xian Science and Technology Press。1992and applicationsM3New York：Springer Verlag，13侯遵泽，杨文采中国重力异常的小波变换与多尺度 1995分析J地球物理学报，1997，40(1)：8595E19 曾昭发，吴燕冈，郝立波磁法技术在浅覆盖区地质HOU Zhun-ze。YANG Wen-caiThe wavelet trans- 填网中的应用R长春：吉林大学，2004formation and multiple scale analysis of Chinese gravi ZENG Zhao-fa，WU Yah-gang，HAO LiboThe ap ty anomaliesJActa Geophysica Sinica，1997，40 plication of magnetic technique in mapping shallow (1)：8595 covered regionRChangchun：Jilin University，143李宗杰小波变换在位场数据处理中的应用I-j石 2004劫劫嘞嘞岛勘嘞勘与物劫物嘞勃岛岛嘞向岛向岛勃勃嘞勘嘞劫嘞岛嘞岛e啤第五届亚太地区地质资源与地质工程国际会议在我校召开10月21日上午9时，由吉林大学、俄罗斯远东国立技术大学、俄罗斯自然科学院普列摩尔斯基分院、 TECHNO-RAPEC(俄罗斯)主办。我校建设工程学院承办的“第五届弧洲太平洋地区地质资源与地质工 程困际会议”在北苑宾馆会议室开幕。我校常务副校长赵继，党委副书记兼副校长韩晓峰，中国科学院院 士林学钰教授，俄罗斯远东国立技术大学教授安德列，本届会议学术委员会主席、我校建设工程学院院长 殷琨，科技处、学科建设办公室、国际合作交流处、地球科学学院、地球探测科学与技术学院、环境与资源学 院、仪器科学与电气工程学院、应用技术学院、水资源与环境研究所等单位负责人参加了开幕式。赵继、安 德烈、殷琨分别致辞。仪式由建设工程学院副院长孙友宏主持。会议为期4天，旨在加快地质资源与地质工程学科的发展，促进世界范围内地质资源与地质工程领域 同行的交流与合作，为环境保护和人类社会的可持续发展建占献策。会议期间，来自俄罗斯、乌克兰、捷 克、也门及国内近百名专家学者将围绕新能源(天然气水合物、油页岩、煤气层、地热等)的勘查评价，远东 及太平洋地区地质资源的开发利用，矿产资源的形成规律与勘探评价，资源开发的环境保护与生态恢复治 理，地下水污染与处理，复杂地层钻进技术，地下信息探测理论与技术，地质灾害勘查、评价、防治与监控新 技术，岩土工程设计与施工及其他地质资源与地质工程的新理论与新技术等问题进行研讨并提交了论文。万方数据基于小波分析的重磁数据求导方法及应用作者：方东红， 曾昭发， 陈家林， FANG Dong-hong， ZENG Zhao-fa， CHEN Jia-lin作者单位：方东红,曾昭发,FANG Dong-hong,ZENG Zhao-fa(吉林大学地球探测科学与技术学院,长春,130026)， 陈家林,CHEN Jia-lin(吉林省区域地质矿产调查所,长春,130022) 刊名：吉林大学学报（地球科学版）英文刊名：JOURNAL OF JILIN UNIVERSITY(EARTH SCIENCE EDITION) 年，卷(期)：2008,38(6)参考文献(19条)1.谢靖 磁源重力异常及其导数的计算方法 1978(03)2.Stavrev P Y Euler deconvolution using differential similarity transformations of gravity or magnetic anomalies外文期刊 19973.曾昭发;吴燕冈;郝立波 磁法技术在浅覆盖区地质填图中的应用 20044.Fishman G S Monte Carlo:concepts,algorithms,and applications 19955.Thomas A;Ridsdill-Smith The application of the wavelet transform to the processing of aeromagnetic data 20006.Torrence C;Compo G A practical gui