地质雷达处理解释的一点体会 中国矿业大学（北京校区）杨峰 地质雷达应用得到广泛的发展和应用。本文主要论述近几年来在地质雷达处理、解释等应用方面的其中一部分做一些探讨，希望得到广大专家的指导。 1 水平预测滤波 从地质雷达采集的信号来看，存在普遍的水平同相轴信号的干扰。这些水平干扰信号对不同深度信号影响效果不同，越深的信号影响越大。这主要是由于高频电磁波在地层传播过程中存在指数形式的衰减和能量扩散等。因而，在相对变化较小的水平干扰信号的作用下，深部反射信号信噪比明显低于浅部信号的信噪比。因此对水平干扰信号的去除，就显得非常重要。如何去除水平信号，其实方法也很多，常用的方法有：（1）背景道去除；（2）窗口滑动平均高通滤波；（3）二维滤波；（4）二维谱的反变化等。不同的方法都有各自的有点，同时也有各自的缺点。不同地区、不同仪器、不同勘探目的、不同采集方法可能都有不同的方法选取。作者在研究去除水平信号过程中，也尝试了不同的数据处理方式，在大量数据试验的基础上，提出水平预测滤波，将信号预测和滤波结合在一起达到去除水平信号的目的。 一、常规不同水平信号去除方法的对比 1．原始信号 2、背景去除 背景去除中背景噪声的计算是对背景道范围进行求均值运算得出的。这种处理方式对由仪器本身或偶合差异引起的噪声具有较好的效果。 采用背景去噪。背景去噪是一种常用的处理方法，尤其对数据量大的剖面运算，速度快，操作简单。但是：在进行背景去噪之前，一定要将剖面上没有的数据（如：天线停滞采集的数据，隧道的蔽塞洞等）删除掉，避免这些信号对综合背景信号的干扰。 3、窗口滑动平均高通滤波 滑动平均高通滤波，其实并不是真正意义上的滤波处理，其原理与背景去除相同，无非该方法的背景噪声是随道数窗口移动，对局部信号的突出有更显著的效果。 4、二维滤波 二维滤波就是利用F－K域中视速度的不同来提取滤波因子，从而达到压制水平信号的目的。 5、二维谱反变化 首先计算出信号的二维谱值，在对二维谱进行编辑，将去除干扰信号的谱值清楚，在通过反变化达到压制干扰信号的目的。由于二维谱运算量较大，这种方法比较适合短剖面的处理。 二、水平预测滤波 水平预测滤波只需要输入一个参数即可，即预测步长。通过不同的预测步长就可以达到水平信号压制程度不同的目的。预测步长越小，其水平信号压制能力越强，否则相反。 2 城市区域雷达勘探高压线缆信号的识别与去除 在城市区域进行地质雷达勘探，无论屏蔽天线还是非屏蔽天线都会受到高压线缆的干扰，如果不对这些干扰信号识别和去除，可能将干扰信号错误地解释为地层信号，甚至整个剖面都不能解释。这里讨论的识别和去除并不是针对所有的信号都有效，这里只是抛砖引玉，希望大家能提出更多的方法，扩大我们的思路。 一、高压电缆信号的识别 当天线从高压电线杆经过时，在雷达剖面上会形成一个清晰的双曲线异常，如果不做分析，可能会将该异常解释为地下金属管线。其实识别的方法很简单，就是利用反射双曲线弧度来求出异常双曲线的速度，根据速度大小，来确定异常的来源。如果所求出速度大于0.2m/ns，
地质雷达学习资料 一．雷达理论基本要点 1.1地质雷达的波组特征 雷达天线发射的是子波而不是单脉冲，子波由几个震荡波形组成，占有一定的时间宽度，反射与折射波依然保持有原来子波的特点，只是幅值上有所变化。这里将雷达子波的周期、持续时间长度和衰减比三个参量作为子波的波阻特征。子波的频率成分与天线的主频相近，持续一个半到两个周期，后续振相略有衰减。例如对于100MHz天线的子波，持续时间可到15-20ns，对于1GHz的天线，持续时间约2ns。子波的波形的确定对于后期处理是非常重要的，它是小波处理的基础。有很多方法可以获得各种频率天线的子波，最简单的方法是利用金属板反射。将一块较大的金属板放置于地面上，发射与接受天线与金属板平行，相距为3个周期的时程，进行数据采集，即可获得子波记录。不同类型的雷达、不同型号的天线，雷达子波的形状是不同的。天线与介质的距离、介质的电导特性对子波的形态和特点也有一定的影响，应根据现场工作条件从记录中分离子波。从下边的记录中也可以辨认出子波的特征。表面反射波、内界面反射波都是近联各州其的衰减波形。对其进行分析可以得到子波的波组特征 为获得雷达探测的结果，需要对雷达记录进行处理与判读，判读是理论与实践相结合的综合分析，需要坚实的理论基础和丰富的实践经验。雷达记录的判读也叫雷达记录的波相识别或波相分析，它是资料解释的基础。在此首先介绍波相分析的基本要点。 1.2雷达波资料解释三要素 要点1：反射波的振幅与方向 从反射系数的菲涅耳（Fresnel）公式中可以看出两点，第一点，界面两侧介质的电磁学性质差异越大，反射波越强。从反射振幅上可以判定两侧介质的性质、属性；。第二点，波从介电常数小进入介电常数大的介质时，即从高速介质进入低速介质，从光疏进入光密介
7地质雷达记录的波相识别 地质雷达反射记录的波形比地震波复杂的多，一方面是由于地质雷达分辨率高记录的信号丰富，另一方面是由于电磁波的干扰因素多，此外还由于雷达发射的子波比较复杂，并非简单的脉冲。因而雷达资料的处理与解释是一项复杂细致的工作。特别是各种地层、目标体、干扰波的识别需要坚实的理论基础和丰富的实践经验。 7.1 地质雷达的波组特征 雷达天线发射的是子波而不是单脉冲，子波由几个震荡波形组成，占有一定的时间宽度，反射与折射波依然保持有原来子波的特点，只是幅值上有所变化。这里将雷达子波的周期、持续时间长度和衰减比三个参量作为子波的波阻特征。子波的频率成分与天线的主频相近，持续一个半到两个周期，后续振相略有衰减。例如对于100MHz天线的子波，持续时间可到15-20ns，对于1GHz的天线，持续时间约2ns。子波的波形的确定对于后期处理是非常重要的，它是小波处理的基础。有很多方法可以获得各种频率天线的子波，最简单的方法是利用金属板反射。将一块较大的金属板放置于地面上，发射与接受天线与金属板平行，相距为3个周期的时程，进行数据采集，即可获得子波记录。不同类型的雷达、不同型号的天线，雷达子波的形状是不同的。天线与介质的距离、介质的电导特性对子波的形态和特点也有一定的影响，应根据现场工作条件从记录中分离子波。从下边的记录中也可以辨认出子波的特征。表面反射波、内界面反射波都是近联各州其的衰减波形。对其进行分析可以得到子波的波组特征 7.2 地质与工程介质结构及反射特征 雷达的探测对象通常是多界面结构，如各类地层、岩性，松散层、风化层等都是多层结构。隧道中的围岩、初衬、二衬等，也是多界面结构。雷达波向介质内传播时，被称为下行波，经反射回表面的波称为上形波。下行波每遇到一个界面就发生一次反射和折射，入射波能量即被分成两部分，一部分经折射继续向下传播，另一部分经反射掉头向上，变成上行波。反射与折射能量的分配与反
GSSI软件RADAN地质雷达资料处理步骤 视图—工具栏（前四个打勾） 视图---状态栏（打勾）状态栏是屏幕显示窗口最下面的那一栏 剖面的线扫描图是每道的波形压缩成一条直线然后用颜色显示出来的。 1.打开软件RADAN，选择文件夹View?Customize（自定义）?Directories（存 放数据的地方）. 2.编辑文件属性，去除只读属性。打开文件File?Open(*.dzt)（原始数据）。 此选项为可选项，一般的.dzt文件不是只读的。 3.扫描信息预编辑：利用图标Edit?Select（选择时避开打的标记， 若选定段包含标记，可把标记避开分几次选）, 选择一段扫描剖面，切除多余扫描信息Cut，或者保存特定扫描剖面 Save。 4.文件测量方向掉转。打开文件，选择File?Save As ，打勾。 5.（针对连测）添加距离信息。测量轮测量直接获取距离概念。连续测量方式 加距离需要三步A) 编辑文件头内的距离信息Edit?File Header, 扫描/ 米[scans/m], 米/标记[m/mark]，B)编辑用户标记（用中文randan5软件时，标记编辑：1）标记类型user 2）编辑模式：添加 3）在线扫描图上需要标记的地方点一下，扫描数自动变化 4）转换—用户标记都转换为 复合标记），C)并利用距离归一化函数进行处理,Process?Distance Norm, ,打勾。（文件前、后没有被标记进去的，归一化后自动去掉。所以标记时一定要从有用的信息开始标，尾部有用信息结束处也要做标记。Usermark用户手标，diatmark测量轮自动打标。） 6.添加里程信息.A)Edit?File Header ?3D option?X start输入里程起点坐标. B)Edit?Edit database?regions ? x origGlobal输入里程起点坐标 ----apply.（此步骤可随时做） 7.水平刻度调整。 Process?Horizontal scale.叠加stacking、抽道 skipping（显示长剖面的整体效果）、加密stretching（波形比例加大，使小的目标体显示清楚）。 8.确定地面反射波信号位置Edit?File Header?position(ns)（把垂直刻度 调成时间显示，然后点地面位置（波尖），状态栏右下角有数值显示，把此数值的负值填到头文件的position里，然后实际深度从调整后的0向下算）。 9.调整信号延时信息，找到地面Process?Correct Position?delta pos (ns). 10.设置和修改介电常数，计算深度信息Edit?File Header?DielConstant。 （测量时的介电常数是虚的，但测出来的剖面图像是真的。通过软件改介电常数，然后软件自动算出来的深度是跟此介电常数相匹配的深度。介电常数求法：1.经验值2、钻孔法，若用钻孔时利用其钻孔深度来修正介电常数-----在头文件里填入不同的介电常数，然后看是否和钻孔真的深度匹配，然后修正介电常数直到和钻孔深度匹配为止。）
几种时频分析方法综述2——希尔伯特—黄变换 夏巨伟 (浙江大学空间结构研究中心) 摘要：希尔伯特—黄变换由经验模态分解（empirical mode decomposition ，简称EMD ）和Hilbert 谱分析两部分组成。经验模态分解方法是一种自适应的、高效的数据分解方法。由于这种分解是以局部时间尺度为基础，因此，它适应于非线性、非平稳过程。通过经验模型分解，任何复杂的数据集都可以被分解为个数有限的、而且常常是为数不多的几个固有模函数(intrinsic mode functions ，简称IMF)的线性叠加。通过分解得到IMF 后，就可以对每一个分量做希尔伯特变换，得到其瞬时频率和幅度。本文详细对Hilbert-Huang Transform 的过程进行了阐述，并用算例分析指出了其优势所在。 关键词：希尔伯特—黄变换；时频分析技术； 1 希尔伯特—黄变换（Hilbert-Huang Transform ） 1.1 希尔伯特变换与瞬时频率（Hilbert Transform and instantaneous frequency ） 对于任意一个时间序列X(t)，它的希尔伯特变换具有如下形式： -1 ()(t)=,-X Y P d t ττπτ ∞∞? 其中，P ——积分的柯西主值； 希尔伯特变换对于任何属于L p 空间中的函数都存立，即上式中X(t)∈L p （— ∞，+∞）。 通过上述定义，X(t)和Y(t)成为一组复共轭对，同时能够构造一个实部和虚部分为X(t)和Y(t)的解析信号(Analytic Signal)Z(t)，Z(t)表示为： ()()(t)=(t)(t)=a ,i t Z X iY t e θ+ 其中， ()()1/222 (t)a =(t)+(t),arctan .X(t)Y t X Y t θ????= ????? 理论上讲有无数种方式去定义虚部，但是希尔伯特变换是唯一能够得到解析 信号结果的方法。 X(t)的Hilbert 变换实质上是将X(t)与函数1/t 在时域上做卷积，这就决定了通过X(t)的Hilbert 变换能够考察其局部特性。得到X(t)的瞬时相位函数后，其瞬时频率为： ()() (t).d w t dt θ= 1.2 经验模态分解与固有模态函数（Empiricalmode decomposition/EMD and Intrinsic mode function/IMF ） 固有模态函数需要满足两个条件：（1）极值与零点的数量必须相等或最多相差一个；（2）由局部极大值包络和局部极小值包络定义的平均包络曲线上任何一点的值为0；
第二讲国内外地质雷达技术发展状况（历史与现状） 探地雷达的历史最早可追溯到20世纪初，1904年，德国人Hulsmeyer首次将电磁波信号应用与地下金属体的探测。1910年Leimback和Lowy以专利形式在1910年的专利，他们用埋设在一组钻孔里的偶极子天线探测地下相对高的导电性质的区域，并正式提出了探地雷达的概念。1926年Hulsenbeck第一个提出应用脉冲技术确定地下结构的思路，指出只要介电常数发生变化就会在交界面会产生电磁波反射，而且该方法易于实现，优于地震方法[1,2]。但由于地下介质具有比空气强得多的电磁衰减特性，加之地下介质情况的多样性，电磁波在地下的传播比空气中复杂的多，使得探地雷达技术和应用受到了很多的限制，初期的探测仅限于对波吸收很弱的冰层厚度（1951，B.O.Steenson，1963，S.Evans）和岩石和煤矿的调查（J.C.Cook）等。随着电子技术的发展，直到70探地雷达技术才重新得到人们的重视，同时美国阿波罗月球表面探测实验的需要，更加速了对探地雷达技术的发展，其发展过程大体可分为三个阶段： 第一阶段，称为试验阶段，从20世纪70年代初期到70年代中期，在此期间美国，日本、加拿大等国都在大力研究，英国、德国也相继发表了论文和研究报告，首家生产和销售商用GPR的公司问世，即Rex Morey和Art Drake成立的美国地球物理测量系统公司（GSSI），日本电器设备大学也研制出小功率的基带脉冲雷达系统。此期间探地雷达的进展主要表现在，人们对地表附近偶极天线的辐射场以及电磁波与各种地质材料相互作用的关系有了深刻的认识，但这些设备的探测精度、地下杂乱回波中目标体的识别、分别率等方面依然存在许多问题。 第二阶段，也称为实用化阶段，从20世纪70年代中后其到80年代，在次期间技术不段发展，美国、日本、加拿大等国相继推出定型的探地雷达系统，在国际市场，主要有美国的地球物理探测设备公司（GSSI）的SIR系统，日本应用地质株式社会（OYO）的YL－R2地质雷达，英国的煤气公司的GP管道公司雷达，在70年代末，加拿大A－Cube公司的Annan和Davis等人于1998年创建了探头及软件公司（SSI），针对SIR系统的局限性以及野外实际探测的具体要求，在系统结构和探测方式上做了重大的改进，大胆采用了微型计算机控制、数字信号处理以及光缆传输高新技术，发展成了EKKO Ground Penetrating Radar 系列产品，简称EKKO GPR系列。瑞典地质公司（SGAB）也生产出RAMAC 钻孔雷达系统，此外，英国ERA公司、SPPSCAN公司，意大利IDS公司、瑞典及丹麦也都在生产和研制各种不同型号的雷达。80年代全数字化的GPR问世，具有划时代的意义，数字化GPR不仅提供了大量数据存储的解决方案，增强了实时和现场数据处理的能力，为数据的深层次后处理带来方便，更重要的是GPR 因此显露出更大的潜力，应用领域得以向纵身拓展。 第三阶段，从上个世纪80年代至今，可称为完善和提高阶段。在此期间，GPR技术突飞猛进，更多的国家开始关注探地雷达技术，出现了很多探地雷达的研究机构，如荷兰的应用科学研究组织和代尔夫大学，法国_德国的Saint-Louis 研究所（ISL），英国的DERA，瑞典的FOA，娜威科技大学和地质研究所，比利时的RMA，南非的开普敦大学，澳大利亚昆士兰大学，美国的林肯实验室和Lawrence Livermore国家实验室以及日本的一些研究机构等等。同时，探地雷达也得到了地球物理和电子工程界的更多关注，对天线的改进、信号的处理、地下目标的成像等方面提出了许多新的见解。GSSI公司在商业上取得了极大的成功，
（一）简单的EMD程序 function imf = emd(x) % Empiricial Mode Decomposition (Hilbert-Huang Transform) % imf = emd(x) % Func : findpeaks x = transpose(x(:));%转置 imf = []; while ~ismonotonic(x) %当x不是单调函数，分解终止条件 x1 = x; sd = Inf;%均值 %直到x1满足IMF条件，得c1 while (sd > 0.1) | ~isimf(x1) %当标准偏差系数sd大于0.1或x1不是固有模态函数时，分量终止条件 s1 = getspline(x1);%上包络线 s2 = -getspline(-x1);%下包络线 x2 = x1-(s1+s2)/2;%此处的x2为文章中的h sd = sum((x1-x2).^2)/sum(x1.^2); x1 = x2; end imf{end+1} = x1; x = x-x1; end imf{end+1} = x; % FUNCTIONS function u = ismonotonic(x) %u=0表示x不是单调函数，u=1表示x为单调的 u1 = length(findpeaks(x))*length(findpeaks(-x)); if u1 > 0, u = 0; else, u = 1; end function u = isimf(x) %u=0表示x不是固有模式函数，u=1表示x是固有模式函数 N = length(x); u1 = sum(x(1:N-1).*x(2:N) < 0); u2 = length(findpeaks(x))+length(findpeaks(-x)); if abs(u1-u2) > 1, u = 0; else, u = 1; end function s = getspline(x) %三次样条函数拟合成元数据包络线 N = length(x); p = findpeaks(x); s = spline([0 p N+1],[0 x(p) 0],1:N);
第34卷第2期福州大学学报(自然科学版)Vol.34No.2 2006年4月Journal of Fuzhou University(Natural Science)Apr.2006 文章编号:1000-2243(2006)02-0260-05希尔伯特-黄变换谱及其在地震信号分析中的应用 陈子雄,吴琛,周瑞忠 (福州大学土木建筑工程学院,福建福州350002) 摘要:介绍了希尔伯特-黄变换(HHT)这一非线性、非平稳信号处理方法,并利用HHT处理了地震工程中 常用的El Centro地震波,得到了该信号的Hilbert谱、边际谱和能量谱,提取了该信号的主要动力特性,并与 该信号的Fourier分析结果进行了对比,显示出HHT这一方法的优越性. 关键词:希尔伯特-黄变换;经验模态分解;固有模态函数;地震信号 中图分类号:TU311.3文献标识码:A Hilbert-Huang transform spectru m and its application in seismic signal analysis CHEN Zi-xiong,W U Chen,ZHOU Rui-zhong (College of Civil Engineering and Architecture,Fuzhou University,Fuzhou,Fujian350002,China) Abstract:HHT is a ne w method to deal with non-linear and non-stationary data.El Centro earth- quake wave is analyzed by HHT.Through the way,Hilbert spectrum,marginal spectrum and energy spec trum are got and dynamic property is extrac ted.The comparison between HHT spectrum and Fourier spec trum is made and the superiority of HHT is demonstrated. Keyw ords:Hilbert-Huang transform;empirical mode decomposition;intrinsic mode function;seismic signal 地震信号具有短时、突变等特点,是一种典型的非平稳随机信号,必须对其进行分析与处理,才可以提取信号的主要特征.传统的Fourie r变换能够表述信号的频率特性,但不提供任何时域信息[1],而小波分析虽然在时域和频域都具有很好的局部化性质,但本质上仍是一种窗口可调的Fourier变换,在小波窗内的信号必须是平稳的,因而没有根本摆脱Fourier分析的局限[2].小波基的选择也是信号分析中的一个重要问题,另外,小波基的有限长会造成信号能量的泄漏,使信号的能量-频率-时间分布很难定量表述. Hilbert-Huang变换(HH T)的信号处理方法被认为是近年来对以Fourier变换为基础对线性和稳态谱分析的一个重大突破[2].它由经验模态分解(E mpirical Mode Decomposition,E MD)方法和Hilbert变换(H T)两部分组成,其核心是E MD分解.该方法采用了固有模态函数(Intrinsic Mode Function,I MF)概念以及将任意信号分解为I MF组成的思想,即E MD法,使得瞬时频率具有实际的物理意义[3].它不受Fourier分析的局限,可依据数据本身的时间尺度特征进行模态分解,分解过程中保留了数据本身的特性,再对各I MF分量进行Hilbert变换,得到信号能量在时间尺度上的分布规律,实现地震动力特性的提取. 1Hilbert-Huang变换 1.1经验模态分解和固有模态函数 经验模态分解(EMD)的目的是通过对非线性非平稳信号的分解获得一系列表征信号特征时间尺度的固有模态函数(I MF),使得各个I MF是窄带信号,可以进行Hilbert分析.首先设定两个条件:1整个时间序列的极大极小值数目与过零点数目相等或最多相差一个;o时间序列的任意点上,由极大值确 收稿日期:2005-07-27 作者简介:陈子雄(1981-),男,硕士研究生;通讯联系人:周瑞忠,教授. 基金项目:教育部博士点专项科研基金资助项目(20040386004)
一般地质雷达数据处理步骤
分界面厚度变化时可用此法，一般不用2）有倾斜地层时可用此法3）使钢筋显示更清楚用此法⑹主要用此法的地方1）测工字钢个数，埋深，形态，间隔2）测空洞3）测钢筋网个数 1.反褶积、一维频率滤波（取默认值。垂直方向上出现一串时（等间隔的多次 波）用此）。Process→Deconvolution；Process→IIR Filter. 2.偏移归位Process→Migration，选择偏移类型kirchhoff，调整曲线形态。 3.希尔伯特变化Process→Hilbert Xform，选phase显示瞬态相位信息。 4.添加地面高程信息，并利用高程归一化函数进行处理。Process→Surface Norm。 5.静态校正Process→Static，mode选择manual手动调整方式。 6.文件拼接。打开Radan软件，选择File→Append files。 7.通道合并，多通道资料对比分析。打开Radan软件，选择File→Combine channels。 8.交互式解释View→Interactive，生成*.lay文件。 步骤1）点2）如果从没解释时就选generate new pick file，如果是在原来的基础上对此文件进行解释就选pick file找到lay文件3）选目标体（如钢筋类的， 解释后可以看出有多少根）：①在剖面上点右键---target options—new target—双击目标体名字----然后在target parameters里改各个要改的参数②在剖面上 点右键---pick options---在pick options里填参数（若拾取工具选block时，在剖面上选一块然后点右键然后加点，）4）选分层①在剖面上点右键----layer options---改layer options里的参数然后确定②在剖面上点右键---pick options---在pick options里填参数（若拾取工具选block时，在剖面上选一块然后点右键然后加点，若当中有空的没有连起来则点右键，插值）5）在剖面上点右键----spreadsheet（表格）6）在剖面上点右键----save changes---current file---保存为lay文件7）用excel打开此lay文件（打开时分割符号选tab键和逗号），打开后去掉头文件然后画图。 速度的选取：在剖面上点右键---ground truth（钻孔）----z（分界面距地面的埋深） 9.绘制地质剖面图.利用电子表格Excel或者Surfer 8软件绘制地质图件。 一：连接文件 File----append files----把每个文件双击------done 二：单个文件宏处理 1）打开文件 2）New macro---保存为宏文件cmf
用希尔伯特黄变换（HHT）求时频谱和边际谱 1.什么是HHT？ HHT就是先将信号进行经验模态分解（EMD分解），然后将分解后的每个IMF分量进行Hilbert变换，得到信号的时频属性的一种时频分析方法。 2.EMD分解的步骤。
EMD分解的流程图如下：
3.实例演示。 给定频率分别为10Hz和35Hz的两个正弦信号相叠加的复合信号，采样频率fs=2048Hz的信号，表达式如下：y=5sin(2*pi*10t)+5*sin(2*pi*35t) (1)为了对比，先用fft对求上述信号的幅频和相频曲线。 代码: function fftfenxi clear;clc; N=2048; %fft默认计算的信号是从0开始的
t=linspace(1,2,N);deta=t(2)-t(1);1/deta x=5*sin(2*pi*10*t)+5*sin(2*pi*35*t); % N1=256;N2=512;w1=0.2*2*pi;w2=0.3*2*pi;w3=0.4*2*pi; % x=(t>=-200&t<=-200+N1*deta).*sin(w1*t)+(t>-200+N1*deta&t<=-200+N2*det a).*sin(w2*t)+(t>-200+N2*deta&t<=200).*sin(w3*t); y = x; m=0:N-1; f=1./(N*deta)*m;%可以查看课本就是这样定义横坐标频率范围的 %下面计算的Y就是x(t)的傅里叶变换数值 %Y=exp(i*4*pi*f).*fft(y)%将计算出来的频谱乘以exp(i*4*pi*f)得到频移后[-2,2]之间的频谱值 Y=fft(y); z=sqrt(Y.*conj(Y)); plot(f(1:100),z(1:100)); title('幅频曲线') xiangwei=angle(Y); figure(2) plot(f,xiangwei) title('相频曲线') figure(3) plot(t,y,'r') %axis([-2,2,0,1.2]) title('原始信号')
地质雷达软件RADAN用户手册美国地球物理测量系统公司 美国劳雷工业公司 2010年10月
RADAN处理软件安装 安装采集软件RADAN66和RADAN5，并且激活采集软件 输入软件序列号serialnumber 输入处理软件产品ID代码:radan 计算获取软件激活码 Windows7系统安装radan5 安装radan程序，找到setup.exe鼠标右键要求以系统管理员身份运行； RADAN软件第一次运行要以系统管理员身份打开。 Windows7系统调整显示效果 选择控制面板->所有控制面板项->显示->更改配色方案->windows经典->高级，对话框如下： 选择颜色 项目->桌面->颜色->设置红绿蓝
资料整理 1打测量，布置网格和测线，数据采集 2数据拷贝与备份： 从地质雷达主机把数据复制在个人电脑上，并利用2种以上存储介质对原始数据进行备份。 3野外记录整理： 整理野外记录本（包括各种参数，利用数码相机或者扫描仪 对原始纪录扫描拍照，并制作成PDF格式文件便于日后随时查看野外现场原始资料），工作照片，收集的各种第三方资料（设计图纸、设计厚度、第三方检测资料），现场钻孔资料（里程桩号、芯样实物和照片、长度）。 利用钻孔资料反算电磁波传播速度或者材料介电常数。 4数据编辑与初步整理RADAN 5资料处理RADAN 6资料解释 7图片制作 8探测报告编写
IGSSI地质雷达探测资料处理流程图数据备份,资料整理,资料处理,资料解释
IIGSSI处理软件功能模块介绍 基本工具 打开数据文件，显示雷达数据剖面。 保存数据文件，保存雷达剖面。 选择数据块，选择目标数据剖面。 剪切数据块，切除多余数据剖面。 保存数据块，单独保存雷达数据剖面。 复制剖面图像至剪贴板，地质雷达剖面制作图片。 编辑数据文件头，输入相关参数：标记间隔、扫描数、介电常数、信号位置； 编辑标记信息、补充漏打的标记、删除多余标记信息。 线扫描显示方式、以灰阶图或者彩色图形式显示雷达剖面。 波形加变面积显示方式、以波列图形式显示雷达剖面。 线扫描+波形显示方式 波形显示方式，以波形方式显示雷达剖面，做频谱分析，显示雷达工作频段。 三维数据体显示 时间切片显示方式 多通道显示方式、多个雷达剖面做比对处理和信息显示。 交互式解释、进行地质解释，绘制地质剖面图、给出钢筋位置。 调整显示参数，设置各种显示参数，保存、调用显示参数文件，设置长标记。 打印雷达剖面 关于RADAN程序版本信息 命令参数按钮
美国地质雷达隧道超前预报工作介绍目前我们国家地下隧道建设工作量大，地质条件复杂，有灰岩地区、花岗岩地区、黄土高原、第四季覆盖等等。 隧道开挖中常常遇到岩溶发育、出现大的空洞，充水或者充泥，有时地下暗河发育；也会遇到构造带，或者岩石破碎，同时地下水发育，这给隧道开挖和建设造成很多困难，同时也给隧道运营造成一定的隐患。因此需要采用一定的手段对这些地质构造和地质灾害进行探测和预报，提前采取措施来排除灾害。 工作任务 为了能够探明隧道开挖面(俗称"掌子面")前方的地质构造，通常采用多种方法进行综合分析、探测、预报。常见的方法有：地质分析，地球物理探测(声波法、直流电法、电磁波方法)，钻孔方法，或者超前导洞等等。采用各种地球物理方法进行探测，分别给出探测结果，综合地质构造情况，进行综合解释，给出掌子面前方的地质构造和可能的地质灾害信息。 探测前提条件 隧道开挖中遇到的地下材料或者介质，主要有石灰岩、花岗岩、大理岩、砂岩、第四季覆盖、沙土、黄土，还有地下水、空洞等等。由于这些材料的物理性质有很多种，比如密度、导电率、介电常数、磁导
率等等。 声波超前预报。由于密度不同、声波传播速度不同，可以采用声波法进行探测，出现了地震波超前预报。 直流电法超前预报。根据导电率的差异采用直流电法，预报掌子面前方材料的导电率差异，尤其是含盐份的地下水表现为良导体、而空气为高阻体； 地质雷达预报。根据导电率、介电常数、磁导率的差异，采用地质雷达高频电磁波方法进行探测，获取掌子面前方材料的介电常数差异信息， 瞬变电磁预报。由于岩石、土壤、水、空气的电磁响应不同，采用瞬变电磁方法探测材料的差异。 目前这4种方法在隧道超前预报中都有使用，尤其是地质雷达超前预报方法得到了普遍使用，利用地质雷达方法在隧道掌子面上进行探测，对隧道开挖超前预报，下面介绍这部分内容。 探测仪器 地质雷达方法通常采用高频电磁波发射法工作，频带范围为几兆赫兹到几千兆赫兹，不同的频率探测深度不同，低频电磁波探测深度较大，因而出现了不同中心频率的天线，商业地质雷达通常采用窄脉冲宽频带电磁波信号工作，一般情况下100兆天线在土壤、破碎的岩石、基岩上探测深度范围从几米到十几米甚至30米左右。 目前隧道开挖地质超前预报距离正好是要求在十几米到30米左
希尔伯特-黄变换（Hilbert-Huang Transform，HHT） 0 前言 传统的数据分析方法都是基于线性和平稳信号的假设，然而对实际系统，无论是自然的还是人为建立的，数据最有可能是非线性、非平稳的。 希尔伯特-黄变换（Hilbert-Huang Transform，HHT）是一种经验数据分析方法，其扩展是自适应性的，所以它可以描述非线性、非平稳过程数据的物理意义。 1 HHT简介[贺礼平.希尔伯特-黄变换在电力谐波分析中的应用研究[D].湖南:中南大学,2009]HHT的发展。 1995年，Norden E.Huang为研究水表面波构思出一种所谓“EMD--HSA”的时间序列分析法，通过这种方法他发现水波的演化不是连续的，而是突变、离散、局部的。 1998年，Norden E.Huang等人提出了经验模态分解方法，并引入了Hilbert谱的概念和Hilbert谱分析的方法，美国国家航空和宇航局（NASA）将这一方法命名为Hilbert-Huang Transform，简称HHT，即希尔伯特-黄变换。 HHT是一种新的分析非线性非平稳信号的时频分析方法，由两部分组成： 第一部分为经验模态分解（Empirical Mode Decomposition，EMD）（the sifting process，筛选过程），它是由Huang提出的,基于一个假设：任何复杂信号都可以分解为有限数目且具有一定物理定义的固有模态函数（Intrinsic Mode Function，IMF；也称作本征模态函数）；EMD方法能根据信号的特点，自适应地将信号分解成从高到低不同频率的一系列IMF；该方法直接从信号本身获取基函数，因此具有自适应性，同时也存在计算量大和模态混叠的缺点。 第二部分为Hilbert谱分析（Hilbert Spectrum Analysis，HSA），利用Hilbert变换求解每一阶IMF 的瞬时频率，从而得到信号的时频表示，即Hilbert谱。 简单说来，HHT处理非平稳信号的基本过程是：首先，利用EMD方法将给定的信号分解为若干IMF，这些IMF是满足一定条件的分量；然后，对每一个IMF进行Hilbert变换，得到相应的Hilbert谱，即将每个IMF表示在联合的时频域中；最后，汇总所有IMF的Hilbert谱就会得到原始信号的时间-频率-能量分布，即Hilbert谱。 在HHT中，为了能把复杂的信号分解为简单的单分量信号的组合，在进行EMD方法时，所获得的IMF 必须满足下列两个条件： 1）在整个信号长度上，一个IMF的极值点和过零点数目必须相等或至多只相差一点。 2）在任意时刻，由极大值点定义的上包络线和由极小值点定义的下包络线的平均值为零，也就是说IMF的上下包络线对称于时间轴。
HHT-希尔伯特·黄变换 1998年，Norden E. Huang等人提出了经验模态分解方法，并引入了Hilbert谱的概念和Hilbert谱分析的方法，美国国家航空和宇航局（NASA）将这一方法命名为Hilbert-Huang Transform，简称HHT，即希尔伯特-黄变换。 HHT主要内容包含两部分，第一部分为经验模态分解（Empirical Mode Decomposition，简称EMD），它是由Huang提出的；第二部分为Hilbert谱分析（Hilbert Spectrum Analysis，简称HAS）。简单说来，HHT处理非平稳信号的基本过程是：首先利用EMD方法将给定的信号分解为若干固有模态函数（以Intrinsic Mode Function或IMF 表示，也称作本征模态函数），这些IMF是满足一定条件的分量；然后，对每一个 IMF进行Hilbert变换，得到相应的Hilbert谱，即将每个IMF表示在联合的时频域中；最后，汇总所有IMF的Hilbert谱就会得到原始信号的Hilbert谱。 与传统的信号或数据处理方法相比，HHT具有如下特点： (1)HHT能分析非线性非平稳信号。 传统的数据处理方法，如傅立叶变换只能处理线性非平稳的信号，小波变换虽然在理 论上能处理非线性非平稳信号，但在实际算法实现中却只能处理线性非平稳信号。历 史上还出现过不少信号处理方法，然而它们不是受线性束缚，就是受平稳性束缚，并 不能完全意义上处理非线性非平稳信号。HHT则不同于这些传统方法，它彻底摆脱了 线性和平稳性束缚，其适用于分析非线性非平稳信号。 (2)HHT具有完全自适应性。 HHT能够自适应产生“基”，即由“筛选”过程产生的IMF。这点不同于傅立叶变换 和小波变换。傅立叶变换的基是三角函数，小波变换的基是满足“可容性条件”的小 波基，小波基也是预先选定的。在实际工程中，如何选择小波基不是一件容易的事， 选择不同的小波基可能产生不同的处理结果。我们也没有理由认为所选的小波基能够 反映被分析数据或信号的特性。 (3)HHT不受Heisenberg测不准原理制约——适合突变信号。 傅立叶变换、短时傅立叶变换、小波变换都受Heisenberg测不准原理制约，即时间窗口与频率窗口的乘积为一个常数。这就意味着如果要提高时间精度就得牺牲频率精度，反之亦然，故不能在时间和频率同时达到很高的精度，这就给信号分析处理带来一定 的不便。而HHT不受Heisenberg测不准原理制约，它可以在时间和频率同时达到很 高的精度，这使它非常适用于分析突变信号。 (4)HHT的瞬时频率是采用求导得到的。
目录 ? 1 本质模态函数(IMF) ? 2 经验模态分解(EMD) ? 3 结论 ? 4 相关条目 ? 5 参考文献 ? 6 外部链接 [编辑]本质模态函数(IMF) 任何一个资料，满足下列两个条件即可称作本质模态函数。 ⒈局部极大值(local maxima)以及局部极小值(local minima)的数目之和必须与零交越点(zero crossing)的数目相等或是最多只能差1，也就是说一个极值后面必需马上接一个零交越点。 ⒉在任何时间点，局部最大值所定义的上包络线(upper envelope)与局部极小值所定义的下包络线，取平均要接近为零。
因此，一个函数若属于IMF，代表其波形局部对称于零平均值。此类函数类似于弦波（sinusoid-like），但是这些类似于弦波的部分其周期与振幅可以不是固定。因为，可以直接使用希尔伯特转换，求得有意义的瞬时频率。 [编辑]经验模态分解(EMD)
EMD算法流程图 建立IMF是为了满足希尔伯特转换对于瞬时频率的限制条件之前置处理，也是一种转换的过程。我们将IMF来做希尔伯特转换可以得到较良好的特性，不幸的是大部分的资料并不是IMF，而是由许多弦波所合成的一个组合。如此一来，希尔伯特转换并不能得到正确的瞬时频率，我们便无法准确的分析资料。为了解决非线性（non-linear）与非稳态（non-stationary）资料在分解成IMF时所遇到的困难，便发展出EMD。 经验模态分解是将讯号分解成IMF的组合。经验模态分解是借着不断重复的筛选程序来逐步找出IMF。 以讯号为例，筛选程序的流程概述如下: 步骤 1 : 找出中的所有局部极大值以及局部极小值，接着利用三次样条 (cubic spline)，分别将局部极大值串连成上包络线与局部极小值串连成下包络线。 步骤 2 : 求出上下包络线之平均，得到均值包络线。 步骤 3 : 原始信号与均值包络线相减，得到第一个分量。 步骤 4 : 检查是否符合IMF的条件。如果不符合，则回到步骤1并且将 当作原始讯号，进行第二次的筛选。亦即 重复筛选次 直到符合IMF的条件，即得到第一个IMF分量，亦即
黄锷院士在《On Holo-Hilbert spectral analysis: a full informational spectral representation for nonlinear and non-stationary data》中提出一种高维全息谱分析理论HHSA(Holo-Hilbert spectral analysis) 要理解HHSA方法，首先要了解希尔伯特变换、经验模态分解(EMD)、与希尔伯特-黄变换(HHT)。 学术背景： 在信号处理与频谱分析的目的是要描述信号的频谱含量在时间上变化，以便能在时间和频谱上同时表示信号的能量或者强度。傅里叶频谱并没有告诉我们哪些频率在什么时候出现。因此傅里叶变换无法表现信号频率成分的时变性，因此学术界先后发展出了短时傅里叶变换、窗口傅里叶变换、小波等手段，近似的求信号某一时刻的瞬时频率。 希尔伯特变换： 希尔伯特变换是以著名数学家大卫·希尔伯特(David Hilbert)来命名。通过希尔伯特变换，使得我们对短信号和复杂信号的瞬时参数的定义及计算成为可能，能够实现真正意义上的瞬时频率的提取，因而希尔伯特变换在信号处理上具有十分重要的地位，使得希尔伯特变换具有广泛的工程应用。 但在进一步的工程应用中，希尔伯特变换具有以下缺陷： (1)希尔伯特变换只能近似应用于窄带信号。但实际应用中，存在 许多非窄带信号，希尔伯特变换对这些信号无能为力。即便是 窄带信号，如果不能完全满足希尔伯特变换条件，也会使结果
发生错误。而实际信号中由于噪声的存在，会使很多原来满足 希尔伯特变换条件的信号无法完全满足； (2)对于任意给定时刻，通过希尔伯特变换运算后的结果只能在 一个频率值，即只能处理任何时刻为单一频率的信号； (3)对于一个非平稳的数据序列，希尔伯特变换得到的结果很大 程度上失去了原有的物理意义。 图1 傅立叶、小波与希尔伯特-黄变换对瞬时频率的分辨率 希尔伯特-黄变换： 针对上述的三个问题，黄锷院士在1998年提出希尔伯特-黄变换(HHT)。其基本思想是：讲一个非稳态、非线性的信号分解为若干个稳态信号，在对分解后的信号进行希尔伯特变换，分别求取对应的瞬时频率。 在这里将非稳态、非线性信号分解为多个稳态信号的算法成为经