【基于MATLAB的地震数字滤波器分析9300字（论文）】

基于MATLAB的地震数字滤波器分析摘要对于地震而言，要降低财产损失和人员伤亡，最有效的方法是利用地震探测仪进行预警，地震探测仪的设计理念就是收集地震信号、转换地震波、过滤地震波、分析地震波和发出预警。其中，过滤地震波是至关重要的一个环节，只要将地震波中的噪声尽可能的过滤掉，留下的有效信号更准确，才更有利于后续的研究和分析。对于滤波器而言，对其性能和功能也提出了越来越高的要求，本课题的研究内容为基于MATLAB的地震数字滤波器。主要设计了Kalman和Wiener两种滤波器，仿真模拟地震信号并加入一个15db和10db的白噪声，比较滤波前后的频率和幅值变化，通过计算信噪比得出Wiener滤波器比Kalman滤波器的降噪效果好。关键词：地震信号;卡尔曼滤波器;维纳滤波器;噪声信号目录TOC\o"1-2"\h\z\u引言 1第一章 绪论 11.1课题研究的内容及意义 11.2国内外研究的现状 1第二章数字滤波器 22.1信号处理的发展和应用 22.2数字信号处理器 32.3数字滤波器 3第三章Wiener滤波器 43.1Wiener滤波器的简介 43.2Wiener滤波器的应用 6第四章Kalman滤波器 74.1Kalman滤波器的基本理论 84.2Kalman滤波器的设计原理 10第五章MATLAB与滤波器的设计 115.1MATLAB的基本应用 115.2使用MATLAB语言进行程序设计 12参考文献 18附录A 19引言我国由于位于世界两大地震带之间，受三大板块的时刻挤压，所以是世界上地震灾害最为严重的国家之一，其地震具有高频率、高强度、浅震源等特点，这一列地震特点导致了地震灾害的突发性、严重灾害性[1]。但随着电子工业的飞速发展，我国的地震探测器技术也在不断增强，但对于地震探测仪而言准确的地震波信号是至关重要的一部分。本次实验主要针对地震信号的数字滤波器进行研究，使用MATLAB仿真软件进行仿真，使用卡尔曼和维纳两种滤波方法设计滤波器，输入模拟地震信号并加入不同频率的白噪声，对比滤波前后幅值和频率波形的变化，计算信噪比，最后通过比较信噪比大小得出维纳滤波器的效果更好。绪论1.1课题研究的内容及意义随着技术的不断发展，数字信号处理在许多领域几乎已经取代模拟信号，人们对它转换的精度的要求也越来越高。电子信息类的工业发展非常迅速，对于滤波器而言，它的功能和性能的要求也越来越多。虽然在以前传统的地震勘探仪器中模拟滤波器使用非常广泛，但是他也存在很多问题，一方面在使用过程中会出来偏移现象，另一方面还会存在器件噪声等等，这些缺点在一定程度上影响了对地震勘探的工作。但是现在数字信号处理技术正在飞速发展，数字滤波器以它稳定的性能、简单的操作、传输方便等特点淘汰了传统的模拟滤波器。1.2国内外研究的现状随着国内经济的不断发展，祖国在建设方面的发展速度也越来越快。但地震对人类的威胁还是不可忽略的。为了将地震引起的伤亡和损失降到最低，地震勘探仪器起着不可忽略的作用。随着科技的进步，地震探测仪也经过了长达二十年的发展，已经实现了数字勘探设备。对于地震勘探仪而言数据采集的数字部分是不可或缺的一个模块，只有将地震信号更好的检测并进行滤波，得到去除掉干扰波的有效波，才能使得地震勘探更加精确。地震信号是一种振动信号，他有两个特点，一个是它的超过120db的动态范围，另一个是高分辨率。地震勘探主要分为三个阶段：数据收集、数据处理和地震材料分析[2]。我们可以通过处理地震数据来获得与地震有关的信息。如何去除信号中的干扰噪声，得到有用的地震波形数据是信号处理过程中的重点、难点[3]。目前，对于地震的勘探而言，设计勘探仪器最主流的方法是利用地震法，地震法的原理是利用检波器来采集到所需要检测的地震波，地震波并不会产生后直接被收集起来，而是在地层中首先要经过传播，其次是折射，最后是反射，然后再返回到地表面，之后再对这些数据进行研究和分析。通过以上介绍可知，当前国际主流的地震信号滤波器方案已经很成熟，性能优秀，广泛应用于万道量级及以下的勘探仪器中，能够满足采集需求[4]。当前地震勘探仪器由负责数据采集的采集站、负责电源供给和数据转发的电源站、负责数据汇总的交叉站和负责记录数据和控制站体的主机构成[5]。对于采集站而言，数字滤波器对于它的设计具有巨大的作用。第二章数字滤波器2.1信号处理的发展和应用在大自然中，信号无处不在，信号的种类有很多，可以分为自然产生的和有意制造出来的，也可以分为有用的和没用的等，无论是哪一种信号，人们通过数字信号处理都可以提取出有用的信号。随着电子工业和信息工业的不断发展，数字信号处理技术已经成为了一个极其重要的问题和技术。数字信号对于模拟信号而言具有以下几点优点。灵活性大。在进行程序的编写和修改时，数字滤波器只需要改变部分程序、代码或者数字就可以实现用户的需求。精度高。数字模电路的元件精度要比模拟电路的高很多。可靠性高。数字系统有且只有“0”、“1”两个电平信号，不易被外界影响而改变。采用大规模集成电路，所以其可靠性也相对较高。易于大规模集成。数字滤波电路的工作主要在截至或饱和状态。并行处理等优点。所以在越来越多的领域中，数字信号发展的十分迅速，并且已经慢慢取代了模拟信号。2.2数字信号处理器数字信号处理就是对信号进行检测和滤波，其算法不仅需要利用数字信号处理器，还需要使用不止一个专用集成电路进行处理。数字信号处理器如同数字信号一般，也具备灵活性高、精确性好等一系列的优点。其在信号处理的过程中一般需要进行三个步骤：信号的过滤与检测、参数的提取和估计、频谱分析等，通过这三个步骤才能使信号变得更加有利于人们的识别。2.3数字滤波器在数字信号处理中，滤波器占有极其重要的地位，数字滤波是语音和图像处理、模式识别、谱分析等应用中的一个基本的处理算法[6]。在各个方向的研究中，人们不单单看重信号处理的快速性，同样对实时性也有很大的要求。在对信号进行处理时，无论是过滤、检测还是预测的过程中，都离不开滤波器的使用，数字滤波器的使用意义无非就是一种离散的时间系统可以把有限精度算法实现。数字滤波器是采用有限精度算法实现离散时间系统，通常应用于修改或者改变频域中信号的属性[7]。数字滤波器具有以下的优点：稳定性高。精度高。设计灵活。实现方便。其中数字滤波器具有稳定性高、精度高、设计灵活、实现方便等许多突出的优点，避免了模拟滤波器所无法克服的电压漂移、温度漂移和噪声等问题，因而随着数字技术的发展，用数字技术实现滤波器的功能越来越受到人们的注意和广泛的应用[8]。应用最多的数字滤波器是线性时不变（LinearTime-Invariant，LTI）滤波器，LTI滤波器通常分为有限脉冲响应（FiniteImpulseResponse，FIR）滤波器和无限脉冲响应（InfiniteImpulseResponse，IIR）滤波器[9]。随机信号和随机过程的存在是非常普遍的。首先，每一个确定性信号在经过采集和测量的过程中总会不可避免地引入一些随机信号，这种随机信号可以被看作是一种随机误差，但是如果不存在随机误差，那么也就无法将确定性随机化了。其次，每一种信号都无法避免随机信号的干扰，所以给这种干扰赋予一个非常形象的称呼——噪声。但是对于干扰和噪声我们又有两种不同的概念，我们可以将干扰看作是一种确定信号也可以看作是噪声，也就是随机信号，它的组成可以分为两个部分，一部分是白噪声一部分是确定性信号。那么区分两种信号通常是判断是否为目标信号，如果是目标信号则为噪声信号，如果不是目标信号，那么就是干扰信号。因此对于滤波器而言，一般是由确定信号伴随这干扰信号同时输入，并且输出时应该尽可能的去除掉干扰信号。维纳滤波和卡尔曼滤波就是解决这类问题的两种方法[10]。噪声可以按照功率谱密度进行划分可以分为以下两种。白噪声（WhiteNoise），他也被称作随机信号。色噪声（ColorNoise）。每一种随机信号都可以看作是由白噪声和确定性信号随机组合的。最优滤波主要解决系统的状态或信号的最优估计问题,即由被噪声污染的观测信号求在某种性能指标和某种意义下状态或信号的最优估值器，也叫最优滤波器[11]。最优的定义是相较而言的，一般情况下对于滤波器选用线性滤波器，因为输入的信号通常是线性函数，对于滤波的定义无非是过滤波形，因为在测量信号时会有一些噪声的存在，噪声会影响对于波形的判断，除了噪声仪器的使用也会导致出现一些误差影响信号的测量结果。通常解决滤波的三种方法分别是：Wiener滤波方法。使用维纳滤波的基本方法可以分为两种一种是利用谱分解的经典维纳滤波，另一种是利用Diophantine方程的现代维纳滤波。Kalman滤波方法。基本方法是利用Riccati方程。现代时间序列分析方法应用现代时间序列分析方法，基于ARMA新息模型和白噪声估值器[12]。本次实验主要使用了前两种滤波方法分别对模拟地震波进行滤波，并选取最适合地震波的滤波方法。第三章Wiener滤波器3.1Wiener滤波器的简介维纳滤波方法是由N.Wiener在20世纪40年代初由于研究火炮控制系统的需要提出来的[13]。它使用的是平稳随机过程谱分析，而信号的平稳随机又成为了它的局限，所以不得不存储所有的历史数据，对于非递推滤波器而言，他不仅存储量非常的巨大，它的计算量也是不容小觑的，所以以上种种缺点导致了他存在不方便使用的问题。但是就算经典维纳滤波由上述种种缺点和问题，但它对于滤波器的设计依旧是非常重要的工具和方法。上世纪八十年代后，又出现了现代维纳滤波方法，他对于经典维纳滤波而言是一种新的突破，它不再局限于平稳随机的信号，而是可以处理非平稳非随机，不仅如此，它还可以处理多维信号，它突破了经典维纳滤波方法的局限。维纳滤波的重要特性有：正交性原理。误差信号统计表征为白噪声[14]。假设，存在一个线性系统，h(n)代表单位脉冲响应，x(n)为输入的一个随机信号，它包含噪声w(n)和有用信号s(n)，则：X(n)=s(n)+w(n)(3-1)则输出为：yn=xn为了使输出信号最大程度接近有用信号，称输出信号为有用信号的估计值s(yn=syx(n)=s(n)=w(n)h(n)图3-1维纳滤波器的系统框图滤波的定义就是用当前的和过去的观测值来估计当前的信号：yn=s(预测，顾名思义就是用过去的观测值来估计当前的或将来的信号：yn=s(平滑或者内插，简单点说就是用之前的观测值来推测过去的信号：yn=sn+上图中的估计信号和想要得到的有用信号是不可能完全相同的，在这里用e(n)来表示真值和估计值之间的误差：en=sn虽然e(n)是随机变量，维纳滤波和卡尔曼滤波的误差准则就是最小均方误差准则：E[e2(n3.2Wiener滤波器的应用设计维纳滤波器需要做几个方面的准备，首先需要估计统计特性，然后再了解清楚观测信号和估计信号直接的关系，以及用到的一些函数，最后才能设计出维纳滤波器，这种设计方法设计出来的滤波器又叫做——后验维纳滤波器。维纳滤波器应用领域十分广泛，例如在生物医学信号处理中非常常见的应用就是有关提取诱发脑电信号的相关研究。研究者在提高信噪比时采用了维纳滤波的滤波方法，前后有多名科学家对滤波方法进行了研究和改进。对于频域应用后验维纳滤波而言，他最重要的地方就是需要使用每一次观测信号分解出来的信号和噪声的谱估计，最后可以通过使用设计的滤波器提高信噪比。在白噪声情况下最有滤波算法本质上是一个利用最小二乘法进行最小方差拟合的过程[15]。考虑真实世界里的仪器，对真实信号的响应(response)不是完美的delta函数,而是由真实信号u(t)与某个已知的响应函数r(t)卷积所得的“被弥散化”(smeared)的信号s(t)：st=−∞∞写成傅立叶变换形式为：S(f)=R(f)U(f)(3-9)另外,仪器所测得的信号c(t)还含有噪音成分n(t)，因此c(t)=s(t)+n(t)。c(t)也称为损坏(corrupted)的信号。在没有n(t)时,我们可以使用去卷积算法,利用傅立叶变换，能很容易地计算得到真实信号：U(f)=S(f)/R(f)(3-10)而此时需要考虑有n(t)时，如何得到u(t)。此时，设存在一个最优滤波器∅(t)orΦ(fUf=C(f)Φ(f)R为了尽量接近真实信号，则使差的模平方最小即可：st=−∞∞将（1）式代入（2）式，得到：−∞∞|Sf对上式，考虑到模的平方，有个小技巧：因为信号与噪音无关(uncorrelated)，所以含有交叉项S(f)*N(f)的积分为0.将S(f)和N(f)分离开，展开模平方，得到：−∞∞1Rf我们要找到一个对测量值有意义的实函数Φ(f),使上式的积分核最小化。因此，将其他量视为常数，对Φ(f)求导并令结果=0,同时考虑到Φf(1−Φf)2'=2可以求得：Φf=|S(f此即最优滤波器，或者Wiener滤波器的公式。这个公式里不含U和R，只需要估计S^2和N^2。当噪音可以忽略时Φf→1,当噪音占主导时，第四章Kalman滤波器维纳滤波理论是四十年代美国科学家Wiener等人发现的，它又有另外一个名称是最佳线性滤波理论。一般来说，维纳滤波在使用过程中存在的最大的难题就是需要使用之前所有的数据，但是这个要求应用在需要随时处理的信号而言是很困难的。为了解决实时性的问题，卡尔曼在六十年代的时候将状态空间模型反映到滤波理论中，并设计出一套递推估计算法，也就是我们所说的卡尔曼滤波理论。基本思想是：采用信号和噪声的状态空间模型，利用前一时刻的估计值和现时刻的观测值来更新对状态变量的估计，求出现时刻的估计值[16]。它的出现不仅弥补了维纳滤波的缺点，而且应用范围也更加广泛了，可以应用在实施运算和计算机运算中。随着电子工业的迅速发展，经典维纳滤波方法已经越来越无法满足实际应用的需求。对于滤波算法的要求最好是递推的，这样方便进行实时计算，然后可以处理多变量非平稳随过程或时变系统滤波问题。卡尔曼滤波是一种时域滤波方法，它解决最优滤波问题的理论基础是状态空间模型和射影，然后再通过计算或求解Riccati方程。卡尔曼滤波的优点是滤波算法可递推，这样方便在计算机上实现，它应用在多个场合中，例如时变系统、非平稳信号和多维信号滤波问题等，现代维纳滤波最大程度解决了经典维纳滤波的局限性和其他缺点，但它并不是完美的，它也需要提前精确了解系统模型和噪声方差。4.1Kalman滤波器的基本理论卡尔曼（Kalman）滤波，是一种非常常见的滤波算法，相对于Wiener滤波而言它具有很多的优点。通过不断的预测和实测来修正自己的估计值，最后达到一个理想的平稳状态[17]。状态空间方法是卡尔曼滤波方法的基本特征和关键技术，状态空间模型是卡尔曼滤波器的设计基础，而状态空间模型的基本概念是系统的状态变量，系统状态变量是描述系统特征的n维列向量，它的取值空间为n维欧氏空间。例如，假设有一个做直线变速运动的点，它的状态通常使用位置、速度和加速度三个要素来表示。如果状态与加速度无关，则可以省略加速度这一要素。状态变量对比信号而言，它的概念更广泛。根据不同的情况有不同的定义，比如说可以把信号看成状态或者状态的分量。状态空间模型一方面是状态空间方法基本模型，另一方面是出发点。状态空间模型又包括两种方程其一是状态方程，其二是观测方程。状态方程又称作状态模型，他通常用来描述状态变化的规律；而观测方程又称作观测模型，它用来定义对状态进行线性观测的方程，需要注意的是观测方程一般情况是含有随机观测噪声。对于线性观测的定义是被观测信号可以是状态本身，也可以是状态的一种或多种分量，也可以是状态分量的线性组合。但观测方程也可以非线性的。一般噪声都是在全频段都有分布的，这是与有用信号处于同一频带的噪声就无法通过FFT算法过滤了[18]。在卡尔曼滤波算法中,首先对某时刻系统的扰动和观测误差(即噪声)展开统计,然后在线性状态空间表示的基础上,通过处理含噪的观测信号求得误差最小时真实信号的估计值,从而排除数据中噪声和干扰信息的影响[19]。卡尔曼滤波适用于所有含有不确定因素的动态系统，可以通过使用一些数学建模对系统接下里的动作作出大致的预测，虽然很难保证系统不会受到外界的干扰，但是卡尔曼滤波在绝大多数情况下都可以提高系统预估的精确度，然后就可以清晰的知道系统状态的转移情况。并且卡尔曼滤波还可以有效的利用多个粗糙数据之间的关系。卡尔曼滤波尤其适用于动态系统，因为它存在以下两个优点：它对于内存的要求很低，因为它仅仅需要保留系统上一个状态的数据，而不需要保留之前很久一段的数据。它的运算速度很快，非常适合解决实时问题和应用于嵌入式系统中。卡尔曼滤波其实是这样一个线性高斯系统应用贝叶斯滤波方程而得到的显式解。hiddenxR+1xRxR-1hiddenxR+1xRxR-1YR+1YRYR-1observedYR+1YRYR-1observed图4-1卡尔曼显示解图形卡尔曼滤波在进行滤波时有两个步骤：第一步是预测，系统根据前一个状态的结果对下一个状态进行判断和估计，第二步是更新，根据目前的结果实时更新系统的状态。对于Kalman而言有五个公式是非常常用的。1、状态先验估计xk|k−1这个公式就是我们研究的变量在这一时刻与上一时刻的具体关系。从这个公式也可以看出卡尔曼滤波只能应用于线性系统中的。2、方差先验估计Pk|k这个公式是用来计算估计值与实际值的方差大小，通常可以利用这个公式评估当前估计的误差有没有问题。3、增益矩阵Kk=Pk卡尔曼滤波通过估计当前观测变量的变化趋势来估计某一时刻变量的值，这个增益矩阵就是指规定的这个变化量的大小。4、状态估计校正xk|这个式子是用来求估计的值。5、方差估计矫正Pk|k=[这个公式是用来修正估计值和实际值的方差。卡尔曼滤波解决的是如何从多个不确定数据中提取相对精确的数据。实验前提是这些数据必须满足高斯分布，而理论前提是第一个高斯斑乘以第二个高斯斑从而得到第三个高斯斑，而这最后一个高斯斑就是提取到相对精确的一个数据范围。卡尔曼滤波器（KalmanFilter,KF）类算法是由美国学者Kalman在20世纪60年代初提出的，此类算法是一种最小方差意义上的最优估计方法[20]。与此同时，卡尔曼滤波器也是一种有关最优化回归数据的处理算法。在很多问题的处理方面，它可以算得上是最优的并且效率也是最高的，毫不夸张的讲它是非常有用的。在过去三十多年的时间段里，它的适用领域不仅仅在日常生活中（人脸识别），还存在于导航领域和控制领域，甚至在军事领域（雷达）都被广泛应用。目前卡尔曼滤波算法在各种复杂情况和环境下的目标跟踪问题上任然是研究和应用热点[21]。4.2Kalman滤波器的设计原理卡尔曼滤波器使用状态观测器。状态观测器是针对可观系统，根据输出y和输入u对系统内部状态u进行观测的结构单元。其构图如下：x=Ax+Buxy=Cx+_+_x=Axxyyxxkk图4-2状态观测器结构图原系统方程为：

x=Ax+Buy=Cx(4-7)采用状态观测器的观测系统方程为：x=Ax+y=Cx为保证观测器的limt→∞x−卡尔曼滤波器就是一种状态观测器，只不过它是随即系统的状态观测器，其结构框图如下：KalmanfilterKalmanfilteryvyyeyvyyeuplantuplantSensornoiseProcessnoiseySensornoiseProcessnoisey图4-3卡尔曼滤波器原理的结构框图在输入u和动态系统plant中间会引入过程噪声w，而在输出y和实际测量yv中间会引入传感器噪声v，而卡尔曼滤波器则是根据u，yv求得测量的最优估计ye。第五章MATLAB与滤波器的设计5.1MATLAB的基本应用5.1.1MATLAB的简介MATLAB是一个高效率的数值计算可视化软件，它的全程叫做“矩形实验室”(MatrixLaboratary)，是美国Mathworks公司开发的一种大型数字计算的软件，他不仅应用范围非常广泛，包括工业、电子、医疗、建筑等领域，例如：信号和图像的处理、检测和分析等，他的优点也十分显著，包括开放性优良、灵活性高、使用简单等。MATLAB在世界范围内也成为了一种非常流行的软件，被越来越多的工作人员、老师和使用。MATLAB语言一方面可以被看成是解释性语言，因为它拥有十分强大的语言编程和交互式计算机的环境，操作者不仅可以利用软件发出一些命令使程序进行操作或数字运算，而且也可以通过使用不同的语言对程序进行编写、运行和调试。另一方面MATLAB语言是一种开放性的语言，用户可以通过实时编写和修改程序以此达到自己的目的。5.1.2MATLAB的特点1、MATLAB在科学计算方面具有完善的工具，所有的变量都是矩阵对象，因为采用矩阵运算而不是循环运算，所以速度快。2、最接近通用语言的科学计算语言。3、支持各种语言的扩展，例如python、c语言等等。4、语法简单，便于操作。5、并行实现较容易。5.2使用MATLAB语言进行程序设计使用模拟地震波，在模拟地震波的基础上加入不同分贝的白噪声，再使用两种滤波器对其滤波，分别对比滤波后和原模拟地震波的时间-幅值和频率-幅值波形的变化，计算出信噪比并进行比较分析，选择效果最好的滤波器。（程序见附录A）Matlab的程序主要分为幅值和频率两个部分，在编写程序时分为以下几个步骤：读取地震波。使用dire函数指定路径。使用filename函数将变量导入工作区，以便后续的使用。使用stream随机流将模拟地震波的数据与时间一一对应画出时间-幅值函数和频率-幅值函数。使用data做出矩阵。对函数做傅里叶变换（fft），生成时间-幅值和频率-幅值图。加入随机白噪声。随机设置白噪声的SNR，单位为db。使用randn(size(data1))函数。利用mean函数求噪声平均值，再用原始噪声减去其平均值，得到噪声信号。计算出信号的的强度后计算出噪声的方差。利用噪声方差的平方根/信号的标准差乘以信号生成噪声信号。利用awgn函数，y=awgn(x,snr,'measured')和y=awgn(x,snr)是相同的，在添加噪声之前测量了x的能量。给正弦波信号加入高斯白噪声生成观测信号。定义误差信号为加入的白噪声（也就相当于观测信号减去地震模拟信号。生成图像。使用卡尔曼滤波算法进行滤波。测量更新预测更新预测误差协方差矩阵测量更新预测更新预测误差协方差矩阵预测误差协方差矩阵预测误差协方差矩阵状态量估计状态量估计预测状态量更新K增益矩阵预测状态量更新K增益矩阵图5-1卡尔曼滤波算法实现流程图设置初始化变量后，计算并赋值转移和测量矩阵、系统噪声矩阵、测量噪声